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ABSTRACT 


This Progress Report covering the period of December 1,1992 to June 
1,1993 presents the development of an analytical solution to the heavy 
ion transport equation in terms of Green's function formalism. The mathe 
matical development results are recasted into a highly efficient compute 
code for space applications. The efficency of this algorithm is accom- 
plished by a nonperturbative technique of extending the Green's function 
over the solution domain. The code may also be applied to accelerator 
boundary conditions to allow code validation in laboratory experiments. 
Results from the isotopic version of the code with 59 isotopes present 
for a single layer target material, for the case of an Iron beam 
projectile at 600 MeV/nucleon in water is presented. A listing of the 
single layer isotopic version of the code is included. 



INTRODUCTION 


Future NASA missions will be limited by exposure to space radiations 
unless adequate shielding is provided to protect men and equipments 
from such radiations. Adequate methods required to estimate the damage 
caused by such radiations behind various shields can be evaluated prior 
to commitment to such missions. 

From the inception of the Langley Research Center heavy ion (HZE) 
shielding program (refs. 1-3), there has been a continued, close relation- 
ship between code development and laboratory experiment (ref. 3). Indeed, 
the current research goal is to provide computationally efficient high 
charge and energy ion (HZE) transport codes which can be validated with 
laboratory experiments and subsquently applied to space engineering design. 
In practice, two streams of code development have prevailed due to the 
strong energy dependence of necessary atomic/molecular cross sections and 
the near singular nature of the laboratory beam boundary conditions (refs. 
4-6). The atomic/molecular cross section dependence is adequately dealt 
with by using the methods of Wilson and Lamkin (ref. 7), allowing effici- 
ent numerical procedures to be developed for space radiations (refs. 
6,8-10). Although these codes could conceivably be applied to the labo- 
ratory validation, methods to control truncation and discretization errors 
would bear little resemblance to the space radiation codes attempting to be 
validated. Clearly, a radical reorientation is required to achieve the 
validation goals of the current NASA space radiation shielding program, 
and such an approach is the main thrust of this research and is briefly 
described below. 



A useful technique in spa.ce radiation shieldinq is the use of the impulse 
response Green's function (refs. 11,12), which satisfies the Boltzman 
equation of the form 
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where Gj m reduces to a monoenergetic unidirectional function at the bound 
ary, Sj(E) iB th® stopping power, is the total cross section, and tfjk 
is the inclusive differential cross section. An arbitrary solution to the 
Boltzman equation within a closed convex region can be written as 
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where f m (E ,n , ,r ) is the incident flux at the boundary (ref. 11). Since 
transport problem is formulated in terms of a single Green's function 
algorithm, the validation of the Green's function in the laboratory 
meets the objective of having a space validated code. Since there is 
hope of a Green's function based on an analytical solution of the 
Boltzmann equation (ref. 13), the resulting evaluation of the shield 
properties should be computationally efficient. 

The first step in this process is to develop an equivalent Green's 
function algorithm in one dimension to match the current capability 
in space radiation transport calculation (refs. 6,14). The algorithm 
is based on the cloised form solution to the one dimensional equation 
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for a monoenergetic beam at the boundary. The probality of validation 
for Ne beams of this algorithm (with multiple scattering corrections) 
has already shown good correlation (refs. 5,15), but improvements in 
the nuclear data base are required for achieving higher correlations 
with experiment. If considerations are restricted to multiple charged 
ions then the right hand side of equation (3) can be further reduced to 
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for which a solution is presented below. 


APPROXIMATE GREENS ' S FUNCTION 

Equation (4) can be simplified by transforming the energy into the 
the residual range as 



and defining new field variables as 

^(i,r j) = Sj(E)tj(x,Bt (6) 
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bo that equation (4) becomes 
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The solution to equation (8) may be written as 
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where zeroth order term of equation (11) is 
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and the first order term of equation (11) is 
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The second order terms of equation (11) are 
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with the condition that G^m( x,r j ,r m) are nonzero f° r 
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The third order terms of equation (11) are 
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and similarly for higher order terms. In the above the g's of n arguments 
are given by 
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In terms of above, the solution to equation (4) may be written as 
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for i=l , the denominator of equation (23) 
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and for i>l, the denominator becomes 
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In equation (22), F m (E) is the integral flux at the boundary, and 

is defined as 
roo 

Fm(E) = J& fm(E , )dB t (26) 

Implementation of equation (22) can now be accomplished independent 
of the character of the boundary values ^ (E) and will give accurate 


results for both space and laboratory applications. 



DISCUSSION OF RESULTS 


Values of collision related fluxes for 3 depths of water target 
for a mono- energetic beam of Iron projectile with 600 MeV/nucleon 
are shown in figures 1 through 3 for both the perturbative and 
nonperturbative Green's function methods for the direct comparison 
of the two methods. The first three collision terms and the sum of 
all collision terms of both theories at a depth of 5 cm of water 
are shown in figure l.A through l.H. The differences in the spectral 
shape is due to the simplification of the attenuation term in the 
nonperturbative theory. The nonperturbation terms represnt the 
average spectrum while perturbation theory retains the spectral shape. 
Figures 2. A through 2.H and 3. A through 3.H are the corresponding 
comparison of the two methods at depths of 10 and 15 cm of water. 

Direct comparsions of figures 1 through 3 shows that the sequence Of 
perturbation terms appear to be converging to a result similar to that of 
nonperturbative result . 

The main advantage of nonperturbative methods are in their 
computational efficencies. The computational time required for the 
nonperturbative code is about 10 minutes on VAX 4000 compared 
to 15, 45, 90 minutes for the 1-st, 2-nd and 3-rd collision terms 

of the perturbation solution. 

Figures 4. A through 4.C show the corresponding differential LET 
spectrum using the method of reference 16. The highest LET peak is due 
to the primary beam and the ion fragments. The successive peaks below 
iron are due to lowdr atomic weight fragments. Such LET spectra can 
be compared to experimental measurments directly. 
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A. 1-st term nonpertubation solution at a depth of 5 cm of water 
for a 600 MeV/nucleon Iron projectile. 
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B. 1-st tdrro perturbation solution at a depth of 5 cm of water 
for a 600 MeV/nucleon Iron projectile. 
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Figure l.E. 3-rd term honperturbation solution at a depth of 5 cm of water 
for a 600 MeV/nucleon Iron projectile. 




r O. 000*16 


O. 00030 



Figure l.F. 3-'rd term perturbation solution at a depth of 5 cm of water 
for a 600 MeV/nucleon Iron projectile. 
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G. All terms nonperturbation solution at a depth of 5 cm of water 
for a 600 MeV/nucleon Iron projectile. 
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Figure l.H. All termd perturbation solution at a depth of 5 cm of water 
for a 600 MeV/nucleon Iron projectile. 
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Figure 2.B. l-st term perturbation solution at a depth of 10 cm of water 
for a 600 MeV/micleon Iron projectile. 



Figure 2.D. 2-nd term perturbation solution at a depth of 10 cm of water 
for a 600 MeV/nucleon Iron projectile. 




Figure 2.E. 3-rd term perturbation solution at a depth of 10 cm of water 
for a 600 MeV/nucleon Iron projectile. 
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Figure 2.F. 3-td terrf* perturbation solution at a depth of 10 cm of water 
for a 600 MeV/nucleon Iron projectile. 
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H. All t^rtnfe perturbation solution at a depth of 10 cm of water 
for a 600 MeV/nucleon Iron projectile. 





Figure 3.B. l-8t term perturbation s 
for a 600 MeV/nucleon Ir 
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Figure 3.H. All terms perturbation solution for a depth of 15 cm of water 
for a 600 MeV/nucleon Iron projectile. 
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Figure 4. 


B. Differential LET spectrum for a depth of 10 cm of water for a 600 
MeV/nucleon Iron projectile. 



nnnrt 


c 

c 

c 


c 


EROGttAM GRfiEN_FUNCTION 


PAfcAMETER( 1 J-59 , IK-1 J*I J , It"U+ 1 , 1ERS-301 ) 

COMMON/T ARGBT/NLAY r X LAY , EM , NAT , ATRG ( 5 ) , ZTRG ( 5 ) , DENSTRG ( 5 ) 


SIKISIISS :1S { ^ lit ; 4Si;?S ; 15 } 1 J ’ 

SiKiUliSS &T i “ i i *« s , . st < 1E « S , . «x , **« s t . r , ^brs , * ^ » 

DATA EO , E IL,BIU, EML,EMU, ACC/IJ*U • ,1K,0.,1K*0.,IK.O.,IR*0*#IL*0./ 


c 


NLAY-1 
NAT- 2 


DN-1. 

ATRG ( 1 ) “ 1 • 

ZTRG( 1 )-l • 

DENSTRG ( 1 )-6.68E22 
ATRG ( 2 ) - 1 6 . 

ZTRG ( 2 ) -8 • 

DENSTRG( 2)=3 . 34E22 


INITIALIZE TS ARRAY WITH A DUMMY CALL, 
ROUTINES ARE ACTIVATED 


SO A( J) ,Z( J) , ANU( J) 


C 


C 


C 

C 

C 

C 


C 

C 


C 


C 

C 


DUMMY— TS ( l . , 1 » 1 ) 

OrEN(UNIT»2, FILE- 
OPEN (UNIT-3, FILE- 


FLUX . DAT ’ , STATUS- ' NEW ' ) 
GLUX . DAT ’ , STATUS- ’ NEW * ) 


X-5. 

BMAX-600. 
PROJECTILE IS 


IRON 26,56 WITH AN INDEX Of 57 IN OUR ISOTOPB TABLE 


1NDEXAP-57 

DO 3 11=1, IJ 
ZT( II ) =Z ( I I ) 
3 AT (II) = A (II) 


koont-o 


DO 1 E-0 . , EMAX , 2 . 
KOUNT-EOUNT+1 


EA=E 4 „ 

IF(EA.LT. 1 .E-3) EA-l.E-3 
ET(KOUNT)=EA 
ST ( KOUNT ) -SMAT ( EA , 1 . , 1 • ) 

RT(KOUNT)=RMAT(EA, 1. , I- ) 
CALL FLUX ( PHI ,EA,X, INDEXAP) 


?°KOUN? ! I ) -PH I ( I ) * ST ( KOUNT ) * ANU ( I ) 
2 CONTINUE 


1 CONTINUE 


IF ( KOUNT. NE.IERS) THEN 
PRINT * , ’ VALUE OF IERS 
STOP 


IN PARAMETER STATEMENT IS INCORRECT' 



c 

c 

c 

c 


c 


END IF 

WRITE ( 2 , * ) i<J , At , zt 

WRITB(2»*)ftT,RT # ST,t- 

WRITE OUT GRAPHICAL OUTPUT (3D), BUT DESCALE FLUX BY LET 
bo 5 k=l,kOUNT 
DO 5 1 = 1 , iNDEXAP 

WRITE ( 3 ,*)I,F(R,I)/(ST(K) *ANU( I ) ) , ET ( K ) 

5 CONTINUE 

CLOSE (UNIT* 2) 

CLOSE(UNIT=3) 

END 


c* ************************************************** 

C SUBROUTINE FLUX (PHI,E,X,IP) 

C PARAMETER ( 1 J=59 ,1K=IJ*IJ,IL=IJ+1 ) 

C DIMENSION B0(IJ),G0(IJ,IJ),Gl(IJ,IJ),G2(IJ,IJ),GG(IJ,IJ) 

DIMENSION BIL(IJ f IJ),EIU(IJ,IJ) , EML( IJ, I J) ,EMU( IJ, I J ) 
DIMENSION B( I J , IJ ) , PHI ( IJ ) 


C 

C 


C 

C 


C 

C 

C 


CALL GREEN(X,E,G0,G1,B,G2,GG,E0,EIL,BIU,EML,EMU) 


DO 1 1=1, I J 
PHI ( I ) =0 . 

1 CONTINUE 

CALL F0(E0( IP),tP,PSI,CPSI,EPSI) 

IFIBO(IP) .GT.500. ) 

-PRINT* , ' E0 ,G0 ,PSI,IP=' ,E0( IP) ,G0(IP ,IP) ,PSI,IP 
PHI ( IP ) =G0 (IP,IP)*PSI 

IF( PHI ( IP) . LT . 1 .E-20) PHI ( IP)=0 . 

ALL BUT NEUTRONS ARE CALCULATED FOR NOW (IE,J>1) 

DO 2 J=2 , IP- 1 
CALL F0 ( EIL( J , IP ) , IP , PSI ,CPSIL, EPSIL) 

CALL F0(EIU(J,IP) , IP, PS I ,CPSIU, EPS1U) 

DNOM=CPSIL-CPSIU 

EAVE= (EIL(J, IP)3-EIU(J, IP) )/2 . 

RAVB=RMAT ( E AVE ,A(IP),Z(IP)) 

IfTdNOM^GT.O. ) EAV={ EPSIL-EPSIU) / (CPSIL-CPSIU) 
IF(EAV.LT.EIL(J,IP) ) BAV=EIL( J , IP ) 
IF(EAV.GT.BIU(J,IP) )EAV=EIU(J,IP) 

IF( (BAV-EIL( J, IP) )*(EIU(J,IP) -EAV) .LT.O. ) 

-PRINT *,' EAV ’ , EAV ,EIL(J,IP) ,EIU(J,IP) 
RAV=RMAT(BAV,A(IP) ,Z(IP) ) 
PHI(J)=PHI(J)+(Gl(J,IP)+B(J,IP)* 

- ( RAV-RAVE ) ) * (CPSIL-CPSIU ) 

C 1F(B(J,IP) • GT . 0 • ) PRINT *,' B= ’ , B( J , IP) , J, IP 

CALL F0(EML( J, IP) , IP , PSI , CPSIL, BPSIL) 

CALL F0 ( EMU (J, IP), IP, PSI ,CPSIU, EPSIU) 

PHI ( J ) =PHI (J)+(G2(J,IP) +GG( J, IP) ) * (CPSIL-CPSIU) 
IF(P!II( J) .LT. l.E-20) PHI ( J ) =0 . 

2 CONTINUE 
C 

RETURN 

END 

C 

c *************************************************** 





SUBttoUf 1MB to ( E * 1 1 , fsi , Cfrsi , tePSl ) 

EJVRAMBtBR( t , 1K»1 3* 1 3 , IL=I J+ 1 ) 

DAtA ZP , AE f BO , SIG, SQ2PI/26 . ,56. ,600. , 1.5,2.5066/ 

ZP,AE,BO,SlG DEFINES THE BEAM 
ZP=CI1ARGE 

EO-MBAN ENERGY (MEV/AMU) 

SIC-STANDARD deviation of energy 

’ ibeAm=1 gAossian dist 

: ibeam=o delta dist 

IBEAM=i 

IF(IBEAM.EQ. 1) GO TO 1 
PSI=0. 

IF(E.GT.EO) THEN 
CPSI=0. 

EPSI-O. 

ELSE 
CPSI-1 . 
epsi=eo 
END IF 
RETURN 

C 1 PSI”SMAT (E f A(IP) r^( IP) ) *TEXP (-( ( EO-B) /SIG) * *2/2 . ) 

EPSI=EO*CPSI+SIG*TEXP ( - ( (E-EOJ/SIG) **2/2 . J/SQ2PI 

RETURN 

END 

C *********************************************** 

c 

FUNCTION ERF(X) 

C 

REAL* 8 ERF , X 
REAL* 8 PERF , PX 

C 

PX= ABS(X) 

ERF=PERF ( PX ) 

I F ( X . LT . 0 . )ERF=-ERF 

RETURN 

END 

C ********************************************** 

C 

FUNCTION PERF(X) 

c 

REAL* 8 Eerf , X ^ fc , 

REAL* 16 QX ,QERF ,T, Al , A2 , A3 

C DATA Al,A2,A3/.3480242D0 f -. 0958798D0 , . 7478556D0/ 

QX=X 

T=1.D0/(1.D0+.47047D0*QX) 

OERF=T* (Al+T* (A2FT*A3 ) ) 

PBRF=1 .DO-QERF* EXP(-QX*QX) 

RETURN 

END 
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c 

c 

c 

c 
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************************************************** 
frUNcTtoN gBERO 1 ( feo , 1 Ap , B , iae , x ) 

FUNCTION gFEROI IS FOR MONO B BEAM AMD VELOCITY CONSERVING INTERACTIONS 

REAL NUP,NUF 

FSl0l=0 . 

GPERO 1=0 . 

IF(B.LE.O.) RETURN 

IF ( Z( IAF ) . LT . 1 . ) RETURN 

FRAGZ=TS ( EO , I AF , I AP ) 

AP=A( IAP ) 

AF=A( IAF ) 

ZP=Z ( IAP ) 

ZF=Z ( IAF ) 

NUP=ZP*ZP/AP 

NUF=ZF*ZF/AF 

R0=RMAT(E0,AP,ZP) 

R J=RMAT ( E , AF , Z F ) 

RJL=NUP* ( RO-X ) /NUF 
RJU=NUP*RO/NUF-X 
B JL=EMAT ( R JL , AF , Z F ) 

E JU=EMAT ( R JU , AF , Z F ) 

ETA=( 2 . *NUP*R0- ( NUP+NUF ) * ( X+RJ ) )/(NUP-NUF) 

XF=(X-RJ-ETA)/2. 

XP=(X*RJ*ETA)/2. 

SIGF=XS( EO ,IAF)-TS(E0,IAF,IAF) 

SIGP=XS(BO, IAP)-TS(E0, IAP, IAP) 

IFtX.EQ.O. JPRINT *,' AF,ZF,SIG . 

IF ( XF. LT. 0. .OR.XP.LT. 0. ) GO TO 98 

PSIOl=FRAGZ*NUF*TEXP ( -SIGF*XF-SIGP*XP ) /ABS( NUP-NUF) 

GPERO 1=PSI0 l/( SMAT ( E , AF , ZF ) /AF ) 

98 CONTINUE 
RETURN 
END 


EQlO.j PRINT *,' AF , ZF, SIGP, SIGF , FRAG= ' , AF , ZF, SIGP, SIGF, FRAGZ 




FUNCTION GPER02(EO,IAP,E,lAF,X) 

Q 

C FUNCTION GPER02 EVALUATES THE SECOND TERM IN PERTURBATION 
C SERIES. THE SOLUTION IS ENERGY INDEPENDENT X SECTIONS 
C AND IGNORES THE MOMENTUM SPREAD CAUSED BY FRAGMENTATION. 

C IE, VELOCITY CONSERVING INTERACTIONS. ALSO A MONOENERGBTIC BEAM. 

C 

REAL NUP,NUJ,NUK 


DATA ZERO/-1 .E-V 
C 

IF(IAF.GE.IAP-l) RETURN 
IF(Z(IAF).LT.l.) RETURN 
AP=A(IAP) 

ZP=Z( IAP) 

NUP=ANU( IAP ) 

SIGP=XS(EO,lAP)-TS(EO, IAP, IAP) 

JAF=IAF 

AFJ=A( IAF) 

ZFJ=Z ( IAF ) 

NUJ=ANU( IAF) 

SIGJ=XS(EO, IAF ) -TS( EO, IAF, IAF) 

C 

DO 99 K=IAF+ 1 , IAP- 1 
KAF=K 




; V 




1 ;t 

i ‘ V 
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AFK-A(KAF) 

ZFK-2(KAF) . 

IF(2FK.Lt.l. ) GO to 98 


NUK-ANU(KAF) 
KS(EO,K 


SIGK-XS ( EO,kAF ) -ts ( EO , KAF , KAF ) 

FRAGK-TS ( EO f K AF , I AP ) 

FRAGJ”TS ( EO, IAF , KAF ) 

XJlMNUKMRMAT ( B l AFK , ZFK ) +X ) -NUP*RO ) / ( NUK-NUJ ) 

X J2=NUP* ( RMAT ( E , AP , ZP ) +X-RO ) / ( NUP-NUJ ) 

IFCNUK.GT.NUP) GO TO 66 
IF(NUJ.GT.NUK) GO TO 77 
X JL=0 . 

IF(XJl.GT.XJL) XJL=XJ 1 
XJU=X 

IF(XJ2 .LT.XJU) XJU=XJ2 
GO TO 80 
66 XJL”0. 

IF(XJL.LT.XJ2) XJL=XJ2 
XJU”X 

IF(XJU.GT.XJ1 ) XJU=XJl 
GO TO 80 
77 XJL=0. 

XJU-X 

IFCXJU.GT. XJl) XJU= XJl 
IF(XJU.GT.XJ2) XJU-XJ2 
80 CONTINUE 

XML=^NUP*RO-NUK*?RMAT ( E , AFK ,ZFK)+X)4( NUK-NUJ ) *X JL) 

XKL”(NUP* ( RMAT( B , AP , ZP ) 4X-RO) - ( NUP-NUJ ) *X JL) 

_ XMu” ( UUP * RO-MUK * ( RMAT ( E , AFK , ZFK ) 4X ) + ( NUK-NUJ ) *X JU ) 

XKU=(NUP* ( RMAT( E, AP, ZP) +X-RO) - (NUP-NUJ) *XJU) 
-/(NUP-NUK) 

FRAGJ”TS ( EO, 1 AF , KAF ) 

IFtXMU . LT . ZERO) GO TO 96 
IF(XKU.LT.ZERO) GO TO 96 
IFCXJU.LT. ZERO) GO TO 96 
IF(XML.LT.ZERO) GO TO 96 
IFCXKL.LT. ZERO) GO TO 96 
IFCXJL.LT. ZERO) GO TO 96 

ARGL=SIGP*XML+S1GK*XKL+SIGJ*XJL 

Anril=SIGP*XMU+SlGK*XKU+SIGJ*XJU 

SLOPE*S IG J+ ( ( NUK-NUJ ) * S IGP- ( NUP-NUJ ) * SIGK ) / ( NUP-NUK ) 
SLOPE=ABS( NUP-NUK)* SLOPE # M 

VALUE” ( TEXP C - ARGL ) -TEXP ( - ARGU ) ) /SLOPE 

GPER02=GPER02+FRAGJ*FRAGK*NUJ*VALUE 

96 CONTINUE 

97 CONTINUE 

98 CONTINUE 

99 CONTINUE 

c 

GPER02=GPER02/( SMAT ( E , AFJ , ZFJ ) /AFJ ) 

RETURN 

END 

************** 4 ************************************ 


c 

c 

c 

c 

c 


FUNCTION GPER03(EP,IAP,B,IAS,X) 

niMrTtnN CPER03 IS AN APPROXIMATE EVALUATION OF THIRD COLLISION TERM 
OF?Se?BR?55SaTIONSeS?e? X-SECTIONS ARE CONSTANT AND THE BEAM IS 
MONOENERGETIC. USE IT AS RECUIRED TIL SOME THING BETTER 
COMES ALONG. PROBABLY NOT TOO BAD. 



non 


c 


HEAL NUf»,NUJ,NUK 


GLX3-0. 

IF(Z(1AS) .LB.O. ) RETURN 
IF(Z( IAF) . LT. 1 . ) RETURN 
ZP=Z(IAP) 

AP*'A( IAP ) 

ZS=Z ( I AS ) 

AS«=Z( IAZ) 

NUP=ANU( IAP) 

NUJ'ANU(IAS) 

GPER03-0. 

IA1=IAS+1 

IA1MAX=IAP-1. 

DO 1 kl=IAl, IA1MAX 
AKl=A(kl ) 

ZKl=Z(Kl) 

SIGRl'OCS(BP,Kl)-TS(EF,Kl,Kl) 

SIGP=XS(EP, IAP)-TS(EP, IAP, IAP) 

SIGS=XS(EP, IAS)-TS(EP, IAS, IAS) 

FRAGR1=TS ( EP , K 1 , IAP ) 

IA2=IAS+1 
IA2MAX=K 1- 1 

DO 3 K2=IA2 , IA2MAX 
AK2=A(K2 ) 

ZK2=Z(K2 ) 

NUK=ANU(K2) 

SIGK2-XS(EP,k2)-TS(EP,K2,K2) 

FRAGK2-TS ( BP , K2 , K 1 ) 

FRAGS=TS(EP,lAS,K2) 

RO=RMAT ( EP , AP , ZP ) 

RLJ=RO-X 

RUJ=(NUP*RO-NUJ*X)/NUK 
EUJ=EMAT ( RUJ , AK2 , ZK2 ) 

ELJ=EMAT(RLJ,AP, ZP) 

IF(NUP.GT.NUK) GO TO 7 
RLJ=NUP * RO/NUK- X 
RUJ=RO-NUJ*X/NUP 
ELJ=EMAT ( RLJ , AR2 , ZK2 ) 

EUJ=EM AT ( RUJ , AP , ZP ) 

IF(NUK.GT.NUP) GO TO 7 
RLJ=RO-X 

RU J=NUP * RO/NUK - X 
BLJ=EMAT ( RLJ , AP , ZP ) 

EUJ=EMAT ( RUJ , AK 2 , Z K 2 ) 

7 CONTINUE 

IF(B.GT.EUJ.OR.E.LT.ELJ) GO TO 3 
GLX3=GLX 3+FRAGS* f RAGK2 * FRAGK I 
-*G3 ( SIGS , SIGK2 , SIGKl , SIGP , X ) /( EUJ-ELJ ) 

3 CONTINUE 
C 

2 CONTINUE 
1 CONTINUE 
C 

GPER03=GLX3 
RETURN 
END 

*************************************************** 

FUNCTION Gl ( A, B, X ) 

C 

BB=B 



IFtA.EO.B) BB=1 . 000000 1*B 

Gl«(TBXi>(-A*X)-TRXP(-BB*X))/(BB-A) 
RBtUNN 
BHD 


C 

C 

C 

c 


******** 




FUNCTION G2 ( A, B,C , X ) 


cc=c 

IF(B.EQ.C) cc-l.oooool*c 
BB-B 

IF(A.EQ.B)BB=1 .0000005*B 
G2=(Gl(A,BB,X)-Gl(A,CC,X) )/(CC-BB) 
RETURN 
END 


C 

C 

C 

C 


********************* 




FUNCTION G3(A,B,C,D,X) 


dd=d 

IF(C.EQ.D) DD=1 . 0000 1 *D 


cc=c 

IF(B.EQ.C) CC=1.000005*C 

G3=(G2(A,B,CC,X)-G2(A,B,DD,X))/(DD-CC) 


RETURN 

END 


C FUNCTION g4(A,B,C,D,E,X) 

C 

Efi—te 

IF(D.BQ.B) EB**1 . 000 1 *E 

G4=(G3(A,B,C,D,X)-G3(A,B,C,EE,X))/(BE-D) 

RETURN 

END 


c 

Q *************************************************** 

c 

FUNCTION 1S(E,J,K) 

C PARAMETER (NE=«7 f I J=59 , NUMl=NE* 1 3 , NUM2=NUM1*I J) 

C dimension tst(Ne, ij,ij),xst(Nb,ij),aa(5),zz(5),dd(5),bt(nb) 

DIMENSION ATST ( I J ) , ZTST ( I J ) 

C COMMON/TARGET/NLA Y , XLAY , DN , NAT , ATRG ( 5 ) , ZTRG(5) ,DBNSTRG(5) 


C 


DATA NOLD/0/ 

DATA ET/25 . , 75 . , 150 . r 300 . , 600 
DATA TST/NUM2*0 . / 

DATA XST/NUM1 *0 • / 


1200. ,2400./ 


NEl-NE 

N=1 

IF(NLAY.NE.NOLD) go to bo 
30 IE=1 

38 IF( 1E.GT.NE1 ) GO TO 35 
ET1*=ET(IE) 

IF(E.LT.ETl) GO TO 39 

IE=IE+1 

GO TO 38 

39 IF(IE.EQ.l) GO TO 35 
TS1=TST( IE, J,K) 



n o 


TS^tS2+(TSl-TS2^*(fe-fef (lfe-i) )/(Et( IB)-BT(lB-l) ) 
GO tO 36 

35 IFIIE.BQ.1) TS=TST(IE,J,K) 

If ( lE.GT.NEl ) TS=TST(NEl, J,K) 

36 RETURN 

ENTRY XS(E,J) 

N=2 

IF(NLAY.NB.NOLD) GO TO 60 
40 IE=1 „ .. 

43 IF(B.LT.BT(IE) .OR.IE.GT.NEl) GO TO 44 

IE-IE+1 

GO TO 43 , . _ 

44 IF(IE.EQ. 1 .OR. 1B.GT.ME1 ) GO TO 45 
XSl=XST(IB, J) 

XS»XS2+|xSl-XS2 ) * (E-ET( IB-1 ) )/(BT( 1B)-ET( IE-1 ) ) 
RETURN 

45 IF(lB.EQ.l) XS=XST(IE,J) 

1F( iB.GT.NEl ) XS=XST(IE-1, J) 

RETURN 

feNTRY ANU(J) 

ANU=ZTST( J)**2/ATST(J) 

1F(ATST( J) .EQ.4 . ) ANU=1 . 000 1 
IF(ZTST( J) .EQ.O. ) ANU= .9999 
RETURN 

ENTRY A{J) 

A=ATST( J) 

RETURN 


'NEWNUC59.DAT' , STATUS 5 * ' OLD’ ) 


ENTRY Z(J) 

Z=ZTST( J) 

RETURN 

80 CONTINUE 
90 OPEN (100, FILE* 

REWIND( 100) 

997 CONTINUE 

READ( 100 , * , END=994 ) IJTST, ATST, ZTST 
9999 FORMAT ( 15/4 ( 8F8 . 2/) /4 ( 8F8 . 2/) ) 

PRINT 9999 , IJTST , ATST , ZTST 
READ ( 100 , * , END=994 ) NN , AA f ZZ , DD 
READ! 100. *» END**994 ) ET,XST,TST 

990 FORMAT( I5/5F10 . 3/5F10 . 1/5E12 . 3/7F10 . l/870( 7E13 . 3/) ) 

PRINT * t ’ TS FOR ' , NN, ' ELEMENTS OF CHARGE ’ , ( ZZ ( I ) , I“1 »NN) 

IF( IJ.NE. I JTST) GO TO 997 
IF(NN.NE.NAT) GO TO 997 

DO 995 JJN=1,NN 

I F ( AA ( J JN ) . NE . ATRG ( J JN ) ) GO TO 997 
IF(ZZ(JJN) .NE.ZTRG(JJN) ) GO TO 997 
IF(DD(JJN) . NE . DENSTRG( JJN) ) GO TO 997 

995 CONTINUE 

PRINT *,' MATERIAL FOR TS ’ , NLAY , ' FOUND oM 4 ^LBAR JIL* ’ 
THERE ARE '.NAT,' ELEMENTS OF CHARGES , ( ZTRG( 1 ) , I“1 , NAT) 

GO TO 996 

" 4 PRINT N *!' MATERIAL FOR TS '.NLAY,' NOT FOUND ON NUCLEAR FILE 
STOP 

996 CONTINUE 

CLOSE( 100, STATUS=’ KEEP’ ) 



nnn 


c 

c 

c 

c 

c 


MoLD-MLAY 8 u 
GO to ( 30, 40)N 
teM) 

SUBROUTIBB GRBB(H*,B,G0,G1,B,G2,GM,B0,EIL,E1U,BML,EHU) 

PARAMETER ,BE-1 .tJ- 5R,RUMl-BEOJ,MUM2-RUMlM J ,IJM-I.J*IJl 

SiUISiiSS l?Ui i :?S‘, c J sS? iit : ! 15 1 J ' 

DIMENSION 

CALL GNFR(X,G0,Gl,B,G2,GM) 

If (B.LT. l.B-10) E-l.E-10 

DO 1 M=1 » I J 

IFCZ(M) .GT.O. ) THEN 

EO(M)=X+RMAT(E,A(M) ,Z(M) ) 

EO(M)='EMAT(EO(M) ,A(M) ,Z(M) ) 

ELSE 
EO(M)-B 
END IF 

DO 1 

eU. j J m -A«Si ' imrmat ( B , M J) ,* 1 J n 

E1U( J iM)*XTANU(J ) *RHAT( E, A( J) • z ( J) )/AHU(M| 

ELSE 

EIL(J,M)=E 

EIU( J,M)=E 
END IF 

rai^iiE^.MJ.AJM, MMjj 

EMU ( J t M ) =EMAT ( E IU ( J » M ) ,A(M) ,Z(M) ) 

ELSE 

EML( J,M)=E 
EMU(J,M)=E 

inEML(J,M) .GT.EMU(J,H) ) THEN 
AAAAB=EML( J pH) 

EML( J,M)=EMU( J,M) 

EMU( J,M)=AAAAB 
END IF 

E1L( J,M)=EML( J,M) 

EIU(J,M)=EMU(J,M) 

S J“SMAT ( E , A ( J ) , Z ( J ) ) J ) 

IF (Z(J) -EQ.O. )SJ=1- 
G0( J,M)~GO( J,M)/SJ 
G1(J,M)-G1( J,M)/SJ 

B(J,M)=B(J,M)/SJ 
G2(J,M)=G2( J,M)/SJ 
GM ( J , M ) =GM ( J , M ) /S J 
1 CONTINUE 
C 

RETURN 

END 

C **************************************** 

C SUBROUTINE GNFR( XA,G0 ,Gl , B ,G2 ,GG) 

C PARAMETER ( I J=59 , I JM=I J* I«J ) 



non 


C 


c 


c 

c 


c 


BIMBNSIOM 

blMBMStoM &( 13 > 13 ) 

X s )tA 

IFtXA.Lt. l.fe-5) X-l.E-5 
CALL GPEHT ( GO , G 1 , G2 , X ) 
CALL GNF(X,GG) 

bO 1 I-lrIJ 


DO 1 J*1 f 1 J 

GG(I, J)=GG( 1, J)-Gl( I , J)-G2( I f J) 

B(I,J)=0. 

B( 1 , 3 ) =ANU 1 1 ) * ANU ( J ) *TS ( 2000 .,I,J)*(G0(J,J) -GO ( I , 
/I f ANUt J1 -ANU( I ) ) *ABS( ANU( 3) -ANU( I ) ) *X) 


1/((ANU(J) 

END IF 
IF (GG(I, J) .LB.O 
IF(I.EQ.J) THEN 


) )*ABS(ANU(J) 
)GG(I,J)=0, 


GG(I, J)=0. 

GO TO 1 
END IF 

DN1=XMANU(J)/ANU(I)-1. ) 

DN2=X*(ANU(J)/ANU( I ) - 1 • ) 
Git I, J)=Gl(I,J)/ABS(DNl) 
G2(I,J)»G2(I, J)/ABS(DN2) 
GG(I,J)~GG(I, J)/ABS(DN2) 
IF( XA . EQ. 0 . ) THEN 
G2(I,J)=0. 

GG(I,J)=0. 


END IF 
1 CONTINUE 


I) 


) 


RETUHN 

END 






C 

C 

C 


C 


C 


C 

C 


SUBROUTINE GNF(X f G) 

PARAMETER ( IJ=59, IL=12, IK=IL*IJ*IJ) 

COMMON/TARGET/NLA Y , X LA Y , DN , NAT , ATRG ( 5 ) , ZTRG ( 5 ) , DENSTRG ( 5 ) 

DIMENSION G( IJ, IJ) ,GT( IL, IJ, IJ) ,G0( IJ, IJ) ,Gl( IJ, IJ) ,G2( IJ, IJ) 
1 , XT( IL) 

DATA XT/0. ,1., 2 . ,4. , 6 . , 8 . ,14. ,20. ,34. ,48. ,62. ,96./ 

DATA NOLD, IPT,GT/0 , - 1 , IK*0 . / 


E0=2000 . 

IF (NLAY.NE.NOLD) GO TO 100 
CONTINUE 

NTAB=IJ**2 „ 

CALL IUNI(IL,IL,XT,NTAB,GT,2,X,G,1PT,IBR) 
IF(IER.NE.0.AND.IER.NE.-4) THEN 

PRINT I UNI IN MODULE GNF FAILED WITH ERROR- ,IER 

STOP 
END IF 

DO 1 1=1, IJ 

DO 1 J=1,IJ 



1 coNtlNOB 


c 

c 


c 

c 

c 

c 


RfttUfcM 

100 CONTINUE 

DO lOl 1 = 1, U 
DO 10 1 J=1,U 

GT(1,I,J)=0. 

IF(I-BQ.J) GT (1,1,J)“1* 

10 1 CONTINUE 

N0LD=NLAY . 

CALL GPBRT(GO,Gl,G2,XT(2)/2. ) 

DO ml i = i»ij 

G0(ll J)=G0( I, J)+Gl ( I,J)+G2(I, J) 

1111 CONTINUE 

DO 102 1=1, IJ 

DO 102 J=1,IJ 
GT(2, I , J)=0. 

GT( 2^1 , J ) =GT ( 2 , 1 , J ) +G0 ( I , K ) *G0( K, J) 
102 CONTINUE 


C 

C 


C 

C 


C 

C 


C 

C 


DO 103 1=1, 1J 

DO 103 J=1,IJ 
GT ( 3 , 1 , J ) “0 . 

GT ( 3? 1 , J ) =GT ( 3 , 1 , J ) +GT ( 2 , I , K ) *GT ( 2 , K , J ) 

103 CONTINUE 

DO 104 1=1, IJ 
DO 104 J=1,IJ 

GT(4,1,J)=0. 

GT(4^1, J)=GT(4,I,J)+GT(3,I,K)*GT(3,K, J) 

104 CONTINUE 

DO 105 1=1, IJ 

DO 105 J=1,IJ 
GT(5,1,J)=0. 

GT(5?I,J)=GT(5,I,J)+GT(3,I,K1*GT(4,K,J) 

105 CONTINUE 
DO 106 1=1, IJ 
DO 106 J=l,lJ 

GT( 6 , I , J) =0 . 

GT(6^I,J)=GT(6,1,J) 4GT(5,I,K)*GT(3,K,J) 
106 CONTINUE 


QWCUm : c , 

OF POOR QUALITY 



DO 107 1*1,1.! 

bo lo7 .1*1,1.! 
gt(7,i,7)=o. 

GT( 7^1, 7,1,7) +GT ( 5 , 1 , K ) *GT( 6 , K , J ) 

107 CONTINUE 

DO 108 1=1,17 

DO 108 7=1,17 
GT(8,I,7)=0. 

DO 108 K=1,I7 _ , _ „ 

GT ( 8 , 1 , J ) =GT (8,1,7) +GT ( 5 , I , K ) *GT( 7 , K ,7 ) 

108 CONTINUE 

DO 109 1=1,17 

DO 109 7=1,17 
GT (9,1,7) =0 . 

DO 109 k=l,17 

GT ( 9 , 1 , 7 ) *GT (9,1,7) +GT ( 7 , 1 , K ) *GT( 8 , K , 7 ) 

109 CONTINUE 

DO 110 1=1,17 

» 

DO 110 7=1,17 
GT( 10 , 1 , 7 ) =0 . 

GT( 10?I^7)=GT( 10,1 ,7) +GT( 7 , I ,K) *GT(9,K,7) 

110 CONTINUE 

"» 

DO 111 1=1,17 
DO 111 7=1,17 

GT( 11, I ,7)=0. 

C 

GT (11, 1^7 ) =GT (11,1,7) +GT ( 7 , I , K ) *GT ( 10,K,7) 

111 CONTINUE 

C DO 112 1=1,17 
C 

DO 112 7=1,17 
GT( 12 , I ,7 ) =0 . 

C 

GT( 12^ I^7)=GT( 12,1,7) +GT ( 9 , I ,K) *GT( 1 1 ,K,7) 

112 CONTINUE 

C . . 

DO 200 1X=1,1L 

C 

DO 200 1=1,17 

c 

DO 200 7=1,17 

1FI1.GT.7) GO TO 200 

IFII.LT.J) GO TO 201 

GT( IX , I ,7 ) =ALOG(GT ( IX , I , 7 ) ) 

GO TO 200 
201 CONTINUE 

1F(TS ( E0 , I ,7 ) . LE . 0 . )GO TO 200 
IF(IX.EQ.l) GT ( I X , 1 , 7 ) =ALOG (TS(E0,1,7)) 


OftS’SlfvAL rAGE SS 

OF POOR QUALITY 



noon non 


tE(IX.GT.l) GT(lX,I, J)=ALOG(GT(IX,I,J)/XT(IX)) 

CONTINUE 

GO TO 2 
END 

,*♦*****♦******************* **************** 

SUBROUTINE GPERT(GGO , GGl ,GG2 , X ) 

PARAMETER! 1*1=59 ) 

DIMENSION GGO(IJ,U) ,GGl(IJ,tJ) f GG2 ( I J, IJ) 

E0=2000 . 

DO 1 

DO 1 J=l,tJ 

IF^I • EqI * GGO ( I , J ) =TEXP ( - (XS(E 0 , I ) -TS(EO , J , J) ) *X) 

1 CONTINUE 

DO 2 1=1 , IJ 
DO 2 J=1,IJ 

-XS( BO , I ) -TS ( EO , I f I ) ,x) 

2 CONTINUE 

DO 3 1=1,13 

DO 3 J=1,1J 
GG2(I,J)=0. 

IF(I.GE.J) GO TO 3 

TO 2 II K j)-GG 2 (I,J| 4 TS(EO,I,E)*TS(EO,E,J)*G 2 (*S(EO,J|-TS(BO,J,Jt, 

-XS(EO^ X)-TS(fio' K*X ) , XS(EO» 1 { -TS(EO, I , 1 ) r X) 

3 CONTINUE 


RETURN 

END 


C ****************************** 
( fuMction BUILD(A,B,X) 




IF ( ABS( A-B ) . LT . 000 1 * A)GO TO 1 
BUILD= ( TEXP( -A*X ) -TEXP ( -B*X ) ) /( B-A) 

RETURN 

1 CONTINUE , 

BUI LD=X* TEXP ( - ( A+B ) *X/2 . ) 

RETURN 

END 

********************************************* 

FUNCTION RMAT(E,A,Z) 

kss,-s.tss.~. ~ 

EXTERNAL ATOE 



non 


c 

c 


SffllwS 1 EHDS»( 8 | ,EROxi(U> ,LU IS | . ENDF1 1 60) ,EHDYH(«0| , 
-feNDXY ( 4 ) ,VAL(6, 1) , WK ( 2790 ) 

Wk SIZE=3*NX*NY+2*NY+NX=3*60* 15+2* 15+60=2790 

COMMON/TARGET/NLAY , XLAY r DM , MAT , ATRG ( 5 ) , ZTRG( 5 ) , DENSTRG( 5 ) 
COMMOM/COM ST /EHDX 1 

DATA lW/0/ , 4 , , i 1 1/ 

DATA iendsw/o, 

DATA IOR/2 , 2/ 

DATA IPl»IP3/-l» » 07 • 08 , . 09 , . 1 » • 2 , • 3 , . 4 , . 5 , 

DATA EM/. 01, .02, .03, .04, .05, . , 0 « i 0 2 0. , 30. , 40. ,50. , 

-60. ,70. ,80. ,92./ 

DATA MOLD, NP/0 , 60/ 

RANGE INTERPOLATIOM/EXTRAPOLATION OF RT(ET,ZT) 

if { Mold . Me . nlAy ) go to iso 

10 CONTINUE 
RMAT»0 . 


IF( E . LE . 0 . J RETURN 
IF(E.GT. .01 JC 


^ ^ l GO TO 20 

?»Li^I(6o'Et,l5,2T,60,B-r,IOB,IPl.Sl.,Z.I.M»T.IBRl 

EAIt-EO ««" ERROR-. «« 

STOP 

Hit I UNI ( 1 5 , 1! 5 . ET , I , END, 1 . 10R.J .1 RVAL, , II P I , IER I 
pr!ot B ;'.' E ?u“ D .h E Sod^e Jeat failed rife error-, ier 

STOP 

END I F 

RMAT=A*TEXP ( RMAT ) / ( Z * Z ) 

RMAT=RMAT*(E/.Ol)**( l.-RVAL) 

RETURN 
20 CONTINUE 

IFCE.LT. 50000.) GO TO 30 
ELP=ALOG( 35000. ) 

SS2*®«!5!!j!» jIT eO.RT.IOJ , IPl.BLP.E.AMATP.IER, 
jRiOT R ;"-'l« SStoJb'rMRI FAILED WIFE ERROR-, IER 
STOP 

cXEl I IBI( 60,ET,15,FT eO.RT.IOR^IRl.EL.E.RHAT.lERI 
PRINT^ * I BI^IN^ MODULE* RMAT^FAl LED HITE ERROR- . IER 
STOP 

RMATP-A.TEXP ( RMAFP)/ < 2* E > 

SSESLS” iSSSJ ® )/ 1 35000 . .REAP ,) /AL0G< 50000 . /35000 . 1 

RMAT=RMAT* ( ( E/50000 . ) * * ( 1 • -RVAL) ) 

RETURN 
30 CONTINUE 



a n o n 


C 

c 


%UD WITH ERROR"' , IER 

STOP 

END if , 

FMAT-A*TEXP ( PM AT ) / ( Z *Z ) 

RETURN 

STOPPING POWER INTERPOLATION/EXTRAPOLATION OF ST(ET,ZT) 

ENTRY SMAT 
H s 2 

IF(NOLD.NE.NLAY) GO TO 150 
40 CONTINUE 
FMAT-0 . 

IF(E . LE . 0 . ) RETURN 
IF(E.GT. . 0 1 J GO TO 50 

«Li L “l(6o!ET,l5,ZT,60,ST.10R,IPl.Et„Z.RMAT.IER| 

SlMD HIT,, ERROR-. 1ER 

STOP 

CALL 1 1 UN I ( 1 5 , 15 , ZT , 1 , ENDX 1 , 1 OR , Z , RVAL, IP3 , IER) 

^lN? R i^‘lUNI D iN E M6K5LE‘ RMAT FAILED WITH ERROR-’, IER 

STOP 
END IF 

RMAT-Z * Z *TEXP ( RMAT ) 

RMAT-RMAT* (E/.01 ) **RVAL 

RETURN 
50 CONTINUE 

IFCE.LT. 50000. ) GO TO 60 
ELP=ALOG( 35000. ) 

CRbL^^BI (EO^EtJ 15, ET, 60, ST, IOR. IPl ,ELP, E ,RMATP, IER) 

PR1MT B *! ,B ' IB^I^MODULE RMAT FAILED WITH ERROR-’. IER 
STOP 

CALl'iBUSO.BT, 15, ET, 60, ST, IOR. IP I .EL, E, RMAT, IER) 

jRIMri^iBI N IN ! MODULB'RMAT FAILBD WITH ERROR-, IER 

STOP 
END IF 

RMATP-Z * Z *TEXP ( RMATP ) 

rmat-rmat”m RMATP-RMAT ) M ALOG1 E) -EL) /( ELP-EL) 

RETURN 
60 CONTINUE 

CALL L IBI ( 60 , ET , 1 5 , ZT , 60 , ST , IOR , IP 1 , EL, Z , RMAT , IER ) 

PRINT R * RE * IBI^I^MODULE RMAT^FAILED WITH ERROR-’, IER 

STOP 
END IF 

RMAT=Z*Z*TEXP(RMAT) 

RETURN 

DISTOPPING POWBR|/0(E| CALCULATION, INTERPOLATION/EXTRAPOLATION 
OF ST ( ET , ZT ) 

ENTRY DSMAT 
N-3 

IF(NOLD.NE.NLAY) GO TO 150 



non 


70 CONTI MOB 
RMAT-0 . 

ifib.lb.o.) MEturm 

lF(B.GT. *01) GO TO 80 

?ALL L SUTP AR ( 60 , 15 , £T , ZT, 60, ST, IENDSW , ENDX 1 , ENDXN , ENDY 1 , BNDYN , 
ENDXY, 1. , 1,EL,Z,IW,VAL,WK, IER) 

PRINT^*^ SUTPAR N IM MODULE RMAT FAILED WITH ERROR-', IER 

STOP 
&MD I F 

CALL lUNI ( 15, 15, ZT, 1 , ENDX 1 , IOR, Z , RVAL, 1P3 , IER) 
IF(IBR.NB.0.AND.IER.NE.-4) THEM 

PRIHT I UNI IN MODULE RMAT FAILED WITH ERROR- , IER 

STOP 
END IF 

RMATl=VAL( 1,1) 

RMAT 1 =Z * Z *TEXP ( RMAT 1 ) 

RMAT2=VAL(2, 1) 

RMAT-RMAT 1 * RMAT 2 / . 0 1 
RMAT-RMAT* ( E/ . 0 1 ) * * ( RVAL- 1 . ) 

RETURN 
80 CONTINUE 

IF(E. LT. 50000. ) GO TO 90 

CALL L SUTPAR ( 60 J 1 5 , ET , ZT , 60 , ST , IENDSW , ENDX 1 , ENDXN , ENDY 1 , ENDYN , 
ENDXY , 1 . , 1 , EL, Z , IW, VAL, WK , IER) 


PRINT * , ' SUTPAR IN MODULE RMAT FAILED WITH ERROR— , IER 

STOP 

END IF 

RMATl— VAL( 1,1) 

RMATl=Z*Z*TEXP( RMATl ) 

RMAT2-VAL12, 1 ) 

RMAT-RMAT 1 * RMAT 2/50000. 

RMAT-RMAT *5 00 00 . /E 
RETURN 
90 COMTINUE 

CALLLSUTPAR (60,l5,ET,ZT,60,ST,I ENDSW , ENDX 1 , ENDXN , ENDY 1 , BNDYN , 
ENDXY, 1 ., 1 ,EL,Z, IW, VAL, WK, IER) 


IF(IER.NE.O) THEN „ , 

PRINT *,' SUTPAR IN MODULE RMAT FAILED WITH ERROR-', IER 

STOP 
END IF 

RMAT1=VAL( 1,1) 

RMAT 1— Z*Z *TEXP ( RMAT 1 ) 

RMAT2-VAL( 2 , 1) 

RMAT-RMAT 1 * RMAT 2 /E 


RETURN 


ENERGY INTERPOLATION/EXTRAPOLATION ET(ZT,RT) 


100 


C 


ENTRY EMAT 

N-4 

I F ( MOLD . NE . NLAY ) GO TO 150 

CONTINUE 

RMAT-0 . 

RL-E 

IF( RL. LE . 0 . ) RETURN 
R— ALOG( Z*Z*RL/A) 


CALL I UN I (60,60, RT ( 1 , 1 ) , 1 ,BT, IOR, R, ROFZ (I),IP3,IER) 

IFC IER . NE . 0 .AND . IER . NE . -4 ) THEN , 

PRINT *,' I UNI IN MODULE RMAT FAILED WITH ERROR-', IER 



STOt> 

6Mb 16 

iio COMTlMUfi 

CKLL IUNI (l5,15,ZT,l, ROFZ , IOR, Z , RMAT , IP3 , IER ) 

If( IER.NE . 0 . AND. IER.NE . - 4 ) THEM 

print i UN i IN Module rmat failed with error= ,ibr 

STOP 
END 1 1 

RMAT=TEXP(RMAT) * nnnn , 

IF ( RMAT. GE. .01. AND. RMAT. LE. 50000. ) RETURN 
IF(RMAT.GT. 50000. ) GO TO 140 

DO 135 1=1,15 
ROFZ ( I ) =RT ( 1,1) 

135 CONTINUE 

CALL I UNI ( 1 5 , 1 5 , ZT , 1 , ROFZ , IOR , Z , RMAT , 1P3 , IER) 

TFf IER.NE.O.AND. IER.NE. -4 ) THEN 

PRINT .IUNI IN MODULE RMAT FAILED WITH ERROR= ,IER 

STOP 

CALL* IUNI (15,15,ZT,1, ENDX 1 , IOR, Z , RVAL, IP3 , IER) 

IFIIER.NE.O. AND. IER.NE. -4) THEN , T - D 

PRINT IUNI IN MODULE RMAT FAILED WITH ERROR= ,IER 

STOP 
END IF 

RMAT=A*TEXP(RMAT)/(Z*Z) 

RMAT=0 . 0 1 * ( RL/RMAT ) * * ( 1 . / ( 1 • -RVAL) ) 

RETURN 
140 CONTINUE 

ELP=ALOG( 35000. ) 

(60?Et! 15 ,ZT, 60.RT , IOR* IPl , 6LP * Z * RMATP , IBR) 

IFIIER.NE.O. AND. IER.NE. -3) THEN 

PRINT IBI IN MODULE EMAT FAILED WITH ERROR® , IER 

STOP 

t?fin 1 1 

CALL IBI ( 60 , ET, 1 5 , ZT , 60 , RT, IOR, IPl , EL, Z , RMAT, IER) 

IFC IER.NE . 0 . AND. IER . ME . -3 ) THEN , 

PRINT IBI IN MODULE EMAT FAILED WITH ERROR® , IER 

STOP 
END IF 

RMATP=A*TEXP ( RMATP ) / ( Z * Z ) 

S™I-SI^((i5"55!*mMP)/(35000..RHRT()/ALOO(50000./35000.1 

RMAT=50000. * (RL/RMAT) 1 ./( 1 .-RVAL) ) 

RETURN 

C 

150 OPEN( 101 , FILE® 'ATOMIC. DAT* , STATUS® 'OLD' ) 

REWIND f 1 0 1 ) 

190 READ ( 101, 160,END=170)DDNN,NN,AA,ZZ,DD,RT,ST,ET 

PRINT 1 6 ''SmAt'fSu55''!nn?' BLEMEMTS OF CHARGE ' f ( ZZ( I ) , 1-1 ,NN) 

160 ^ATlho ^5/?nS.3/5Fi0.1/5El2.3/100(5El3.4/,/180(5El3.4/)/ 

-18(5B13.4/) ) 

I F ( DDNN . NE . DN ) GO TO 190 
IF(NN.NE.NAT) GO TO 190 

C 

DO 100 1=1, MM 
IF(I.GT.NAT) GO TO 190 
IF(AA(I) .NE.ATRG(I) ) GO TO 190 
IF(ZZ( I ) .NB.ZTRG( I ) ) GO TO 190 
IF(DD(I) .ME.DEMSTRG(I) ) GO TO 190 



noon non 


180 CONTINUE 

PRINT MATERIAL FOR EMAT ' , NLAY , ' FOUND ON ATONIC FILE ' 

l,' three Are * , nat , * elements of charges ',(ztrg(i),i=i,nat) 

GO TO 200 
170 CONTINUE 

print All Atomic data failed rmat test ' 

CALL MGAUSS( 1 . R-6 , EN( 1 ) , 6 , ANS , ATOE , F , 15 ) 

DO 2 10 1-1,15 

V=ZT( I ) 

W=2.*V 

RT(1,I)=ANS(I) 

DO 210 J=1,NP 
ST(J,I)=SIONA(EN(J) ,W,V) 

210 CONTINUE 

DO 220 1=2, NP 

CALL MGAUSS ( EN ( I - 1 ) , EN ( I ) , 6 , ANS , ATOE , F , 1 5 ) 

DO 220 J=1,15 

RT(I, J)=RT(I-1,J)+ANS(J) 

220 CONTINUE 

DO 240 1=1 , NP 
ET ( I ) =ALOG ( EN ( I ) ) 

DO 240 J=1 , 15 
V=ZT( J) 

W«2 • *V 

RT(i,J)=ALOG(V*V*RT(l,J)/W) 

ST( I , J ) =ALOG( ST( I , J ) /( V*V) ) 

240 CONTINUE 



, NLAY, ’ PLACED ON NE1 
OF CHARGES ',(ZTRG( 


NOLD=NLAY 


DO 250 1=1, 15 

DSl=(ST(2,I)-ST(l,I) )/(ET(2)-ET(l)) 

DS2=( ST (3, 1)- ST (2,1) )/( ET( 3 ) -ET( 2 ) ) . , , 

ENDXl(I)=DSl+(DS2-DSl)*(ET( 1)-ET(2) ) /( ET( 3 ) -ET( 1 ) ) 
250 CONTINUE 


C 

PRINT 260 

260 FORMAT ( ' RMAT IS INITIALIZED') 
GO TO ( 10,40,70, 100) ,N 
END 




♦FUNCTION RTIS( E, A, Z )♦♦♦♦♦♦♦♦**♦** #*♦*♦**♦♦ 


FUNCTION RTIS( E , A, Z ) 

FUNCTION RTIS GENERATES THE RANGE, ENERGY, AND LET ARRAYS FOR 
ARBITRARY IONS IN TARGET. E IS ION ENERGY IN MBV/AMU 


C 


EXTERNAL ATOE 

DIMENSION IOR(2) , IP 1 ( 2 ) ,ROFZ( 15) 



non 


c 

c 

c 


btMfeMStoM teT(60) # RT(60,l5) f ST(60,l5),ZT(15) r ANS(15),F(15) 

bffiliSS *S&d 4**I ‘ i 1 1 i* . i«Di» < 1 5 , . B«D* U 6 0 , , ®«D*« C «0 > . 

-ENDXY(4 ) ,VAL(6* 1 ) ,WK(2790) 

WK S1ZE=3*NX*NY+2*NY+NX=3* 60* 15+2*15+60=2790 

COMMON/TARGBT/NLAY , XLAY , DM , NAT , ATRG ( 5 ) , ZTRG ( 5 ) , DENSTRG ( 5 ) 
COMMON/CONST/ENDX 1 

DATA IW/0/ . . t/ 

DATA IENDSW/0 , 1 , 1 , 1 , 1 r 1 r 1 p 
DATA IOR/2 , 2/ 

: iSsi: :Sfc:T»i*. . 

-60. ,70. ,80. ,92./ 

DATA NOLD,NP/0, 60/ 

RANGE INTERPOLATION/EXTRAPOLATION OF RT(ET,ZT) 


N=1 

IF(NOLD.EQ.O) GO TO 150 
10 CONTINUE 
RTI S=0 . 

1FIE.LE.0.) RETURN 
I F ( E . GT ..01) GO TO 20 

IF ( IER.NE - 0 . AND. IER . NE . 3 ) THEN 

' tM ftirtnflT t? RTT<? FA 


IBI IN MODULE ~RT IS FAILED WITH ERROR* 


IER 


PRINT * 

STOP 

CALL* I UN I ( 1 5 , 1 5 , ZT , 1 , ENDX 1 , 10R , Z , RVAL , IP3 , IER ) 
IF(IER.NE.0.AND.IER.NE.-4) THEN httm errors ' IER 

PRINT *,' IUNI IN MODULE RTIS FAILED WITH ERROR 

STOP 

END IF ,, 

RTIS=A*TEXP(RTIS)/(Z*Z) 

RTIS=RTIS*(E/.01)**( l.-RVAL) 

RETURN 

20 CONTINUE 

IF(E. LT. 50000. ) GO TO 30 
ELP=ALOG( 35000. ) 

ck^S!Sni!l5.«,60,W,IOR,IPl.El.P,*.RT.SP.I«l.l 
prS^^’Ir^h'motulr'rI.s FAILED RITE ERROR-. IER 
STOP 

CALL'lBKSO.ET.lS.ET.eO.RT.IOR.IPl.BL.Z, RTIS, IER) 

primt***^ * ibi^*i)** module” rt is** failed WITH error-. IER 


STOP 

RTISP=A*TSXP(RTISP ) /( Z*Z ) 

, / < 3S»00. .RTIS, ,/ALOO, 50000 
RTIS=RTI S* ( (E/50000. ) * * ( 1 • RVAL) ) 


/35000 . ) 


return 

30 CONTINUE 
EL=ALOG(E) 



noon 


c 

c 

c 


CALL iBl(60,Ef,l5,ZT,60,RT,IOR,IPl,BL,Z,RTIS,IER) 

imER.NB.O.AMDilBR.NE.-J) THEM 

PRI^IT 101 IM MODULE RTIS FAILED WITH ERROR”* ,IBR 

STOP 
END IF 

RTIS=A*TEXP( RTIS ) / ( Z*Z ) 

RETURN 

STOPPING POWER INTERPOLATION/EXTRAPOLATION OF ST(ET,ZT) 

ENTRY STIS 

N=2 , rn 

IF(MOLD.EQ.O) GO TO 150 
40 CONTINUE 
RTIS=0. 

IFCB.LB.O.) RETURN 
I F ( E . GT . .01) GO TO 50 

CALL IBI(60,ET,l5,ZT,60,ST, 10R, IPl , EL, Z ,RTIS, IER) 

WITH ERROR-, IER 

STOP 

CALL IUNI ( 15, 15 ,ZT, 1 ,ENDXl , IOR,Z,RVAL, IP3, IER) 

PRINT^^’lUNI^N^^^ WITH ERROR= ' , IER 

STOP 
END IF 

RTIS=Z*Z*TEXP(RTIS) 

RTIS=RTIS* ( E/ . 0 1 ) * *RVAL 
RETURN 
50 CONTINUE 

IF(E. LT. 50000. ) GO TO 60 
ELP=ALOG( 35000 . ) 

CALL^BI ( 60 , ET ! 1 5 , ZT , 60 , ST , IOR , IP 1 , ELP , Z , RTISP , IER ) 

PRIWT^*^ IBI^IH ^MODULE~RTIS^ FAILED HITH ERROR- , IER 
STOP 

CA^IBIteOEET, 1 5, ZT, 60, ST, IOR, IPl, EL, Z, RTIS, IER) 

PRINT^*^ * IBI^IN^MODULE RTI S^FAILED WITH ERROR- , IER 

STOP 
END IF 

RTISP=Z*Z*TEXP( RTISP) 

RTIS=Z*Z*TEXP(RTIS) . 

RTIS=RTIS+( RTISP -RTIS) * ( ALOG(B) -EL) /(ELP EL) 

RETURN 
60 CONTINUE 

CALL L IBI E 60,ET, 1 5, ZT, 60, ST, IOR, IPl, EL, Z, RTIS, IER) 

lF(IER.NE.0.AND.IER.NE.-3) then 

PRINT IBI IN MODULE RTIS FAILED WITH ERROR= ,IER 

STOP 
END IF 

RTIS=Z*Z*TEXP(RTIS) 

RETURN 

D( STOPPING POWER ) /D( E ) CALCULATION, INTERPOLATION/EXTRAPOLATION 
OF ST(ET , ZT) 

ENTRY DSTIS 
N=3 

IF(NOLD.EQ.O) GO TO 150 
70 CONTINUE 



non 


RTIS®0 . 
If(B*Lfe 
if(B.GT 


return 
< . oi 


^ I ) GO TO 8o 

CALL SUTPAR ( 60 , 15 , ET, ZT, 60, St , IENDSW, ENDXl , ENDXN, ENDY 1 , ENDYN, 
ENDXY , 1 . , 1 , EL, Z , IW, VAL, WK, IER) 

PRINT SUTPAR IN MobULE RTIS FAILED WITH ERROR®', IER 

STOP 

HMD If 

CALL IUNI ( 15, 15,ZT, 1, ENDXl , IOR, Z , RVAL, IP3 , IER ) 
IF(IER.NE.0.AND.IER.NB.-4) THEN 

PRINT I UN I IN MODULE RTIS FAILED WITH ERROR® , IER 

STOP 

END If 

RTISl=VAL( 1,1) 

RTISl=Z*Z*TEXP(RTISl ) 

RTIS2=VAL( 2,1) 

RTIS®RTISl*RTIS2/.0l 
RTIS=RTIS* ( E/ . 0 1 ) * * ( RVAL- 1 . ) 

RETURN 
80 CONTINUE 

IF(E. LT. 50000. ) GO TO 90 

C AL^SUTP AR ( 60 1 1 5 , ET , ZT , 6 0 , ST , IENDSW , ENDX 1 , ENDXN , ENDY 1 , ENDYN , 
ENDXY , 1 . , 1 ,EL, Z , IW, VAL, WK , IER ) 

PRINT SUTPAR IN MODULE RTIS FAILED WITH ERROR® , IER 

STOP 

END IF 

RTISl=VAL( 1,1) 

RTISl®Z*Z*TEXP(RTISl ) 

RTIS2=VAL( 2,1) 

RTIS=RTISl*RTIS2/50000 . 

RTIS=RTIS*50000 ./E 
RETURN 
90 CONTINUE 

C AL^SUTP AR ( 6 0 , 1 5 , ET , Z T , 60 , ST , I BNDSW , ENDX 1 , ENDXN , ENDY 1 , ENDYN , 
ENDXY, 1 1 ,EL,Z, IW,VAL,WK, IER) 
iFflfeR 0) THEN 

PRINT i,’ SUTPAR IN MODULE RTIS FAILED WITH ERROR® ' , IER 

STOP 
END IF 

RTISl=VAL( 1, 1) 

RTlSl=Z*Z*TEXP(RTISl ) 

RTIS2=VAL( 2,1) 

RTIS=RTISl*RTIS2/E 

RETURN 


ENERGY INTERPOLATION/EXTRAPOLATION ET( ZT,RT) 


C 


ENTRY ETIS 
N=4 

IF(NOLD.EQ.O) GO TO 150 
100 CONTINUE 
RTIS=0 . 

RL=E 

IF( RL. LE . 0 . ) RETURN 
R=ALOG( Z*Z*RL/A) 


CALL IUNI(60,60,RT(1,I) , 1 , ET, IOR,R,ROFZ( I ) , IP3, IER) 
IFIIER.NE.O.AND. IER.NE.-4 ) THEN 

PRINT I UN I IN MODULE RTIS FAILED WITH ERROR®', IER 

STOP 



ftttb tt 

no cottttMue 

CALL 1UNI(15, l5,ZT, l,ROFZ,IOR,Z,ATIS,lt>3,IER) 

luSl D iM E MODULE RTIS FAILED WITH ERROR 5 * ' , IER 

STOP 
END It - 

5r ( RT1 S*GB R !o 1 ! AND. RTIS. LE.50000.) RETURN 
IF(RTIS. GT. 50000. ) GO TO 140 

DO 135 1-1,15 
ROFZ( I )=RT( 1,1) 

135 CONTINUE 

CHL IUNItlS,lS,ET,l,ROFE,IOR,E,RTlS,IP3,IER| 

fRIN?"i”' IUNI D iN E H6™LE RTIS TAILED HITE ERROR- , IER 
STOP 

CALL 1 IUNI ( 1 5 , 1 5 , ZT , 1 , ENDX 1 , IOR, Z , RVAL, IP3 f IER ) 

rtis failed with brror=, ' ibr 

STOP 
BHD IF 

imS=A*TEXP(RTlS)/(Z*Z) 

RTIS=0 . 01* {Rl>/RT1S ) ** ( 1 • /( 1 • -RVAL) ) 

RETURN 
140 CONTINUE 

ELP=ALOG( 35000 . ) 

CA^M^t!.5.TT,60.RT,IOR,IPl.ELT.E.RTISE.IER, 

«« ERROR-. IER 

STOP 

CALL 1 IBl(60,ET, 15,ZT,60,RT, IOR, IP 1 , EL, Z ,RTIS, IER) 

PRINT R *? E *1BI N IN* MODULE ETIS FAILED WITH ERROR**', IER 
STOP 

END 1F , „ w , 

RTISP**A*TEXP ( RT I SP ) / ( Z * Z ) 

RVAL^ALOGlIsOOOo! .RTISP)/( . * BT I S > |/AIXX5( 50000./35000. 

RTIS*=50000.*(RL/RTIS)**( l./( l.-RVAL) ) 

RETURN 

150 OPEN( 101, FILE- 'ATOMIC. DAT’ , STATUS=’OLD* ) 

190 REAt) R (* 101^160, END= 170) DDNN , NN , AA , ZZ , DD t RT , ST , ET 

1*™1 ELEMENTS OF CHARGE ',<ZZ(I),I= 

160 FORMAT( ,10^3 .15/5,^10 . 3/5F10 . l/5E12.3/IEO(5B13.T/)/lEO,5E13 

- 18 ( 5El 3 . 4/ ) ) 

I F ( DDNN . NE . DN ) GO TO 190 
IF(NN.NE.NAT) GO TO 190 

C DO 180 1=1, NN 

IF(I.GT.NAT) GO TO 190 
IFIAA(I).NE.ATRG(I)) GO TO 190 
IF(ZZ(I) .NE.ZTRG(I) ) GO TO 190 
IF(DD(I) .NE.DENSTRG(I) ) GO TO 190 
180 CONTINUE 


1 ,NN) 
.4/)/ 



nnn o n n 


MtiNt *i» MAfftfliAt FoR RTiS SMLAY,' FoUNd ON Atomic File ' 
if'THERB ARB SNAT,' ELEMENTS OF CHARGES ' , <ZTRG( I ) , 1 = 1 ,NAT) 

GO TO 200 

170 printT' All Atomic data tailed etis test ' 

CALL MGAUSS( 1 . B-6 , EN( 1 ) , 6 , ANS, ATOE, F, 15 ) 

DO 210 1=1 , 15 
V=ZT(I) 

W=2.*V 

ET(1,I)=ANS(I) 

DO 210 J=1,NP 

ST ( J , 1 ) =S IONA ( EN ( «J ) , W , V ) 

210 COMTINUE 

CALL^MGAUSS(BM(I-1) ,EN(I) , 6, AMS, ATOE, F, 15) 

bO 220 J=1 , 15 

RT ( I , J ) =ET ( I - 1 , J ) + AMS ( J ) 

220 COMTIMUE 

I 

DO 240 1=1, MP 
BT( I )=ALOG( EN( I ) ) 

I 

DO 240 J=1 ,15 
V=ZT( J) 

W=2.*V 

RT(I,J)=ALOG(V*V*RT( I, J)/W) 

ST(I,J)-ALOG(ST(I,J)/(V*V) ) 

240 COMTINUE 

OPEM 1 1 0 3 , F I LE= ' MEWRT I S . DAT ' , STATUS- ' MEW ' ) 

WRITE ( 103, 160 )DN, NAT, ATRG, ZTRG,DENSTRG,RT, ST,ET 

PRINT MATERIAL FOR RTIS ’ ,NLAY, ' PLACED ON NBWRT IS r 

1 ’ THERE ARE ' ,NAT, ' ELEMENTS OF CHARGES ,(ZTRG(I),I 1 , NAT ) 
CLOSE ( 103, STATUS-' KEEP' ) 

200 CONTINUE 

CLOSE ( 101 , STATUS- 'KEEP' ) 

C 

NOLD-NLAY 

C 

DO 250 1 = 1, 15 „ . , 

DSl=(ST(2,I)-ST(l,I) )/(ET(2)-ET(l)) 

ENDXlf I j =DSl+f DS2-DS S M ET( 1^-ET ( 2 j ) / ( ET( 3 ) -ET( 1 ) ) 

250 CONTINUE 
C 

PRINT 260 

260 FORMAT ( ' RTlS IS INITIALIZED') 

GO TO ( 10,40,70, 100) ,N 
END 

********** *** + ***O*SUBR0UTINE ATOE ( E , F )***♦♦** ********** ******** ** 
SUBROUTINE ATOE (E,F) 

SUBROUTINE ATOE IS THE INTEGRAND OF RANGE FORMULA 
DIMENSION ZT ( 1 5 ) , F ( 15 ) 

DATA ZT/1 . , 2 . , 3 .5, 5. , 7 . , 10. , 14 . ,22. , 30. ,40. ,50. , 

-60. ,70. ,80. ,92./ 



bo Id 

V-Jt(t) 

W-2**V 

Fm«w/$ioNAtB,w,v) 

10 CONTINUE 

c 

return 

BHD 

C»******************* FUNCTION SIONA( B, A, Z )******** ****************** 
C 

FUNCTION SIOHA( fi f A, Z ) 

Q 

c e IS IN UNITS OF MEV/AMU 
c 

I F ( 2 . EQ . 1 . ) THEM 
CALL ATOPP ( E , S ) 

SIOMA=S 

CALL ATOPN(E, SN,A, Z) 

SIONA=>SIONA+SN 

return 

ELSE 

CALL ATOPA(E,S f Z) . _ 

S IONA^S* (ZEFF(E,A,Z)/ZEFF(E,4 . ,2. ) ) * *2 
CALL ATOPN(E,SN,A,Z) 

SIONA=SIONA+SN 

RETURN 

END IF 

END 

C* ******************** FUNCTION ZEFF(E, A, Z) ************************** 

C 

FUNCTION ZEFF(E f A, Z) 

C FUNCTION ZEFF CALCULATES ION EFFECTIVE CHARGE 
c 

IF(E.LT.l.B-lO) E^l.E-10 
ARG=1 . /( ( 1 • +E/938 . 3 ) * *2 ) 

IFCARG-EQ. 1. ) THEN 

ZEFF“Z*( 1-TEXP( -$ . 77 *SQRT( E) /Z** ( 2 • /3 . ) ) ) 

ELSE 

BBTA=SQRT( l.-ARG) 

ZBFF=Z * ( 1-TEXP( -125. *BETA/Z** ( 2 . /3 . ) ) ) 

END IF 
RETURN 
END 

£*♦********♦*********** SUBROUT I NE ATOPP (T,S)**************‘******* 

C 

SUBROUTINE ATOPP (T,S) 

C SUBROUTINE ATOPP CALCULATES PROTON STOPPING POWER 
c 

DIMENSION A(12, 5), STP(5) 

C COMMON/TARGET/NLAY , XLAY , DN , NAT , ATRG ( 5 ) , ZTRG ( 5 ) , DENSTRG ( 5 ) 

C 

DATA NOLD/O/ 


IF(NOLD.NE.NlAY) GO TO 100 
CONVERT T IN MEV TO B IN KBV 


C 

C 


10 E=1000.*T 
ENERGY 0 <E <*= 10 KEV 



non non non n non 


lR(g*Lft»l0. ) fHEN 

DO 50 1=1, NAt 
STP ( 1 ) =A (1,1)* SQRT ( E ) 
50 CONTINUE 

ENERGY 10 <B < 1000 REV 


ELSE 1F(E.LT. 1000. ) THEN 

DO 70 1=1, NAT 
SL=A( 2,1 ) *E** . 45 

SI1=(A(3, I)/E)*ALOG( 1 . +A( 4 , 1 ) /E+A( 5 , 1 ) * 
STP( !)=!./( l./SLFl ./SH) 


70 CONTINUE 

ENERGY 1000<=E<=100000 KEV 
ELSE 

B2=( V/C) **2 IN ENERGY TERM 


E) 


B2=l.-l./( 1.+E/938300. )**2 
CALCULATE ASYMPTOTIC DENSITY EFFECT 


C 


C 


C 


C 

C 


C 


DBL=DENSEFF(T) 

DO 30 1=1, NAT 
SUM=0 . 

DO 20 <3=1,5 

L=J-1 

K=J+7 

SUM=SUM+ A ( R , 1 ) * ( ALOG ( E ) * * L ) 
20 CONTINUE 


SUM=SUM* ( 1 . -TEXP( -2 . * 1 1 . B5/E) **2 ) ) M . 

STP( I ) = ( A( 6 , 1 ) /B2 ) * ( ALOG( A( 7,I)*B2/(1-B2)) -B2 -SUM-DEL/2 . ) 

CONTINUE 


END IF 
S=0. 

DO 90 1=1, NAT 
S=S+STP( I ) *DENSTRG( I ) 
90 CONTINUE 


S=S* 1 . E-2 1 
RETURN 


C 

100 CONTINUE 
C 

DO 110 1-1,12 

DO 110 J=l,NAl 

A(I, J)=ACOEF( 1 ,ZTRG( J) ) 

110 CONTINUE 
C 

NOLD=NLAY 
GO TO 10 
END 

C***************** FUNCTION DENSEFF(B) ***************** 

C 

FUNCTION DENSEFF ( E ) 



non n non n nnnnnn 


SUBROUTINE DENSfeKK CALCULATES *HE 'ASYMPTOTIC' DENSITY EFFECT 
OF THE STOPPING POWER (CORRECTION TO BETHB-BLOCH FORMULA) 

E IS THE RiNBTIC ENERGY OF PROJECTILE IN MEV/AMU 

COMMOM/TARGET/MLA Y , X LAY , DN , NAT , ATRG ( 5 ) , ZTRG ( 5 ) , DENSTRG ( 5 ) 


DATA NOLD/O/ 

RYDBERG CONSTANT RY AND ELECTRON REST MASS EMASS ARE IN BV 

DATA RY, EMASS, EMASS, AO, PI/ 13. 6, . 5 1 1E6 , 938 . 3 , 5 . 2928E-9 , 3 . 14159265/ 


IF(NOLD.NE .NLAY ) GO TO 5 
3 CONTINUE 

WP=4.*RY*SQRT(PI*RHOE*(AO**3) ) 

C=-2 . *ALOG( EMEAN/WP ) - 1 . 

ATA=P/MC (KINETIC MOMENTUM/REST MOMENTUM) 


ATA=SQRT(E*E+2 . *PMASS*E ) /PMASS 
DEL=ALOG ( AT A* * 2 ) +C 
IF(DEL.LT.O. ) DEL=0 . 

DENSEFF=DEL 
RETURN 
5 CONTINUE 

Q 

C CALCULATE IONIZATION POTENTIAL DUE TO IONIZATION OF ALL 
C TARGET MATERIALS 

C REMEAN=ALOG OF ALL IONIZATION POTENTIALS 
C RHOB=E LECTRON DENSITY OF ALL IONS 

C 

RBMEAN=0 . 

RHOE=0 . 

C 

DO 2 J=1 , NAT 

RHOE=RHOE+ZTRG( J) *DENSTRG( J) ^ 

REMEAN=REME AN+ZTRG ( J ) * DENSTRG ( J ) *ALOG( ( 2 . * EMASS )/ACOEF( 7 , ZTRG( J) ) ) 

2 CONTINUE 
C 

BMBAN=TEXP ( REMEAN/RHOE ) 

RHOE=DN * RHOE 
NOLD=NLAY 
GO TO 3 
END 
C 

C******»************** SUBROUT I NB ATOPA (T,S)******************** 

C 

SUfiROUTlNB ATOPA(T,S) 

C 

C SUBROUTINE ATOPA CALCULATES ALPHA STOPPING POWER 
C DIMENSION A(9,5),B(9,2),C(9,5), STP ( 5 ) 

C COMMON/T ARGET/NLA Y , X LA Y , DN , NAT , ATRG ( 5 ) , ZTRG ( 5 ) , DENSTRG ( 5 ) 

C 

DATA NOLD/O/ 


C 

C 

C 


C 


ARRAY B IS BASED ON LOW ENERGY H r O SOLID 


DATA B/. 9961, .4126,6.92,8.831,2.582,2.371, .5462,-. 07932,- .006853, 

-1.766, .5261,37.11,15.24,2.804,3.782, .3734, -.1011, -.007874/ 


I F ( NOLD . NE . NLAY ) GO TO 110 
10 CONTINUE 



n nnn n nnn non 


C 


DO 30 1=1,9 


DO 20 J=1,NAT 
A(I,J)=C(I, J) 

IF(ZTRG( J) -EQ. 1 . ) A( I , J ) =B( I , 1 ) 
IF(ZTRG(J) .EQ.0. ) A(I,J)=B(1,2) 
20 CONTINUE 

30 CONTINUE 

CONVERT T IN MEV TO E IN KEV 

E=1000.*T*4 . 

ENERGY E< 1 KEV 

1 F ( E . LT . 1 . ) THEN 

DO 60 1=1, NAT 
STP(I)=A( 1 , 1 ) *E* *A( 2 , 1 ) 

60 CONTINUE 


ENERGY 1 <=E<=10000 KEV 

ELSE IF(E. LB. 10000.) THEN 
DO 80 1=1, NAT 

SHIG=A(3| I)/(E/1000. j *ALOG( 1 . +A( 4 , 1 ) /( E/1000 . ) +A( 5 , 1 ) *E/1000 . ) 
STP( I )=SLOW*SHIG/ ( SLOW+SHIG) 

80 CONTINUE 

ENERGY E>=10000 KEV 


C 


C 


C 


ELSE 

FB=ALOG( 1000./E) 


STP( I ) =TEXP^ A( 6,1) +A( 7 , 1 ) *FE+A( 0 , 1 ) *FE* *2+A( 9,I)*FB**3) 
10 CONTINUE 


END IP 

S=0. 

DO 100 1=1, NAT 
S=S+STP( I ) *DENSTRG( I ) 
100 CONTINUE 


C 

S=S*1 .E-2i 

ACO=l . -( E-2 .E4 ) /2 . E4 
IFIE.LT. 2. E4) RETURN 
IF(E.GT. 4 .B4 ) ACO=0 . 
BCO=l . -ACO 
CALL ATOPP ( T , STPB ) 
S=ACO*S+4 . *BCO*STPB 
RETURN 
110 CONTINUE 
C 

DO 120 1=1,9 


DO 120 J=1,NAT 
C( I , J ) =CCOEF ( I , ZTRG( J ) ) 
120 CONTINUE 


C 


NOLD=NLAY 



noon non o o non non non non ooo o noo o noon 


GO TO 10 
END 


C 

C*****i4At***nm***** * * SUBROUTINE ATOPN( T, S, A1 ,Zl)*************** 

C 

SUBROUTINE ATOPN ( T , S , A 1 , Z 1 ) 

SUBROUTINE ATOPM CALCULATES NUCLEAR STOPPING POWER OF 
ARBITARY IONS 

DIMENSION SN(5),A(5),Z(5),STP(5) 

COMMON/TARGET/NLA Y , X LA Y , DN , NAT , ATRG ( 5 ) , ZTRG ( 5 ) , DENSTRG ( 5 ) 

CONVERT T IN MEV TO E IN KEV 

E=l000 . *T*A1 

DO 40 1=1, NAT 
A( I )=ATRG( I ) 

Z ( I ) =ZTRG( I ) 

EPRM IS REDUCED ION ENERGY 

EPRM=32.53*A( I ) *E/( Z 1*Z( I ) * ( Al4A( I ) ) * ( Z 1** ( 2 . /3 . )+Z ( I ) ** ( 2 ./3 . ) ) 

REDUCED ION ENERGY EPRM <.01 

If (EPRM. LT. 0.0 1 ) THEN 
SN( I )»1 . 593* EPRM** . 5 

REDUCED ION ENERGY . 0 1 <=EPRM<=10 

ELSE If (EPRM. LE. 10. ) THEN 

SNf I)=1.7*(EPRM**i5)*( ALOG ( EPRM+TEXP (l.))/(1.46.B* BPRM+3 . 4 * EPRM 
-**1.5) ) 

REDUCED ION ENERGY EPRM>10 
ELSE 

SN(I)=ALOG( .47*EPRM)/(2.*EPRM) 

END IF 
40 CONTINUE 

CONVERT FROM LSS REDUCED STOPPING UNIT TO UNITS OF EV/10**15 ATOMS/CM**2 


S=0. 

DO 60 1=1, NAT 

STP( I ) =SN( I)*(8.462*Zl*Al*Z(I) )/( ( Al+A( I ) ) *SQRT( Z 1** ( 2 . /3 . )+Z(I) 
-**(2./3. ) ) ) 

S=S+STP( I ) *DENSTRG( I ) 

60 CONTINUE 

S=S*1 .E-21 

RETURN 

END 

********************** FUNCTION ACOEF ( I r Z)************ + ******* + ***** 
FUNCTION ACOEF ( I, Z) 

FUNCTION ACOEF CONTAINS ZIEGLER STOPPING POWER COEFFICIENTS 
FOR PROTONS 



DIMENSION A(12,36),ZB(5),AA(12,5) 
bATA ((A(l,J),l“i,l2),j=l,14)/1.2 


-.02972, -.0006243, 3. 469 ,3. 907, 5725. ,1461.,. 008829, 
--9.449, 3.635,-. 5001, . 02961 ,-. 000642 1 , 3 . 519, 3. 963, 
-.007782, .01326,3650. ,-9.809, 3. 763, -.5 164, .0305,-. 
-3. 535, 6288., 1372. , .007361, .01 377, 3453., -10. 17, 3.8 

-.03139,- .0006779, 3 .553, 4 .004 ,6205. ,555.1, .008763, 
-3297. ,-10. 53, 4. 019, -.549, .03229, -.0006957/ 


03506 . - *UUU /D j 9 $ o • 1 t ' 9 f 1 r J • u r • uvu 

-3154., -11. 95, 4. 537, -.6 169, .036 13, -.0007759, 5. 554, 6. 3, 6496., 
-110., .009689, .01632, 3097., -12. 34 ,4. 68, -.6358, .03721, -.00079 
-5.323,6.012,7611. .292.5, .006447, .01683,3024. ,-12.72,4.023, 
—.6547, .03828, -.0008203, 5. 874, 6. 656, 7395., 117. 5,. 007684, 

-.01734,3006. ,-13. 11, 4. 965, -.6735, .03935, -.0008425, 5. 611, 

-6.335,8046. ,365.2, .006244, .01785, 2928., -13. 4, 5.083, -.6906, 

-. 04042, -. 0008675, 6. 4 11, 7. 25, 8262. ,220. , .006007, .01836, 2855. 
— 13. 69, 5. 2, -.7076, . 04 15 ,-. 0000925/ 

DATA ZB/40.,47.,72.,82.,92./ 

DATA AA/6. 734, 7. 603, 9 120. ,405.2, .005765, .02039, 

-2704. ,-14. 59, 5. 463, -.7333, .04242, -.0008990, 
-5.623,6.354,7160; ,337.6, .01394,-02396, 

-2193. .-18. 39, 6. 86, -.92 11, .05346, -.001 14, 


C 



canono non 


J=2 

If(2<Lfe.3M tHBN 
AcoEf"=A( i ) 

RETURN 

ELSE lfr(Z.LT.ZB( 1) ) THEN 

ACOEF=A( I , 36 ) +( AA( I, 1) - A (1,36) )*(Z-ZB(1) )/( ZB( 1 ) -36 . ) 

RETURN 

ELSE 

C 

DO 3 11=1,5 
IIZ=II 

IF(J.LT.ZB(1I)) GO TO 4 

3 CONTINUE 
C 

4 CONTINUE 

ACOEF=AA(l,ltZ-l) + (AA(I,IIZ)-AA(I,IIZ-l) )* 

-(Z-ZB( 1 IZ-1 ) )/(ZB( I IZ ) -ZB( IIZ-1) ) 

RETURN 
END IF 
END 

****************** ** function ccoef( i,zj *************************** 

FUNCTION CCOEF ( I , Z ) 

FUNctlON CCOEF CONTAINS ZIEGLER STOPPING POWER COEFFICIENTS 
FOR ALPHA PARTICLES 

FOR A DISCUSSION SEE HEALTH PHYS. VOL. 46, P.1101,1984 

DIMENSION C( 9 , 37 ) ,CC (9,6) , ZB( 6 ) 

C 

C ARRAY C IS BASED ON THE FOLLOWING LOW ENERGY GAS/SOLID CHOICES 
C G=GAS,S=SOLID 

C ll-G, HE-G, LI-S, BE-S, B-S, C-G , N-G, O-G, F-G, NE-G, NA-G,MG-S, AL-S, 

C SI-S,P-S,S-G,CL-G,AR-G,K-S,CA-S,SC-S,TI-S,V-S,CR-S,MN-S,FB-S, 

C CO-S,NI-S,CU-S f ZN-S ,GA-S ,GE-S, AS-S, SE-S,BR-G,KR-G,RB-S 

C ARRAY ZB IS BASED ON 
C ZR-S,MO-S, AG-S, AU-S,PB-S,U-S 

C 

DATA ((C(I,J) ,1=1,9) , J=1 , 18)/ 

-.39, .63,4. 17,85.55, 19.55,2.371, .5462, -.07932, -.006853, 

-.58, .59,6. 3, 130. ,44 .07,2.809, .4847, -.08756, -.007281, 

-1.42, .49, 12.25,32. ,9. 161,3.095, .44 34, -.09259, -.007459, 

-2.206, .51, 15.32, .25,8.995,3.28, .4 188, -.09564, -.007604, 

-3.691, .4128, 10.48,50.72,9. ,3.426, .4, -.09796, -.0077 15, 

-3.47, .4485,22.37,36.41,7.993,3.588, .3921, -.09935, -.007804, 

-2. , .548,29.82, 1 8 . 1 1 , 4 . 37 , 3 . 759, . 4094 , - . 09646, - . 00766 1 , 

-2.717, .4858,32.88,25.88,4.336,3.782, .3734,-. 1011, -.007874, 
-2.616, .4708,41.2,28.07,2.458,3.816, .3504,-. 1046, -.008074, 
-2.303, .4861,37.01,37.96,5.092,3.863, .3342,-. 1072, -.008231, 
-13.03, .2685,35.65,44 . 18 , 9 . 175 , 3 . 898 , . 3191 , - . 1086, - .00827 1 , 
-4.3, .47,34 .3,3.3, 12.74,3.961, .314,-. 1091, -.008297, 

-2.5, .625,45.7, .1,4.359,4.024, .3113,-. 1093, -.008306, 

-2.1, .65,49.34,1.788,4. 133,4.077, .3074, -.1089, -.008219, 
-1.729, .6562,53.41,2.405,3.845,4.124, .3023, -.1094, -.00824, 
-3.116, .5908,53.71,5.632,4.536,4.164, .2964, -.1101, -.008267, 
-3.365, .571,63.67,6. 182,2.969,4.21, .2936, -.1103, -.00827, 
-2.291, .6284,73.88,4.470,2.066,4.261, .2994,-. 1085, -.00814 5/ 
DATA ( (C(I,J) ,1=1,9) ,J=19,37)/8.554, .3817, 

-83.61, 11.84, 1.875,4.3, .2903,-. 1103, -.008259, 

-6.297, .4622,65.39, 10. 14,5.036,4 .344, .2897,-. 1102, -.00824 5, 
-5.307, .4918,61.74, 12.4,6.665,4.327, .2707,-. 1127, -.00837, 
-4.71, .5087,65.28,8.806,5.948,4.34, .26 18, -.11 38, -.0084 2, 
-6.151,-4524,83., 18 . 31 , 2 . 71 , 4 . 361 , . 2559, - . 1 145, - .00844 7 , 
-6.57, .4 322,84 .76, 15 . 33 , 2 . 779 , 4 . 34 9 , . 24 , - . 1 166 , - . 00855 , 
-5.738, .4492,84 .61, 14. 18,3. 101,4 .362, .2327,-. 1174, -.008580, 



nan 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


-6 »0l3,. 4707,86.66,16. 65, 3. 211,4 *375, .2253, -.1185, -.008648, 

-4*32» *4947,76. 14, io« 86. 5 *44 1.4 *362, .2069, - .1214, • 

-3 114*. 62 36, 76. 67, 7. 62, 6. 38 5, 4. 355, . 18,-. 1255, -.009045, 

-.’.,* ..it ■»« m 7 62 7 • 502 , 4 . 389 , » 1 806 , - • 1253 , — . 009028 , 

"3* l ll'* Wit'll' 67. 7.62:8:514 J^O?!. 1759, -.1258, -.009054, 

-5.746^ ‘4662! 79. 24^ 1. 185, 7. 993, 4. 4 19, .1694, -.1267, -.009094, 

‘ra^-5Sa:ia:5:5”5rj:U^J:S55;:13U;i:135S:=:S855S5: 

I! 6^:7 H8.1, l!47,. 9686, 4. 436,. 1443, -. 1299, -.009229, 

1**13 7377 , 1*7. 9. 1-466, 1.016,** 478, .1600,-. 1262,- 

i.jji:*. 37i, *.*•*.. i5i7. -.iJ7»,-.o«»«7t/ 

SK£ ocJIS: 6*^ :*1T i -SB'* iUS* • 1 ial*“ * , 5o5*i6*^** t ' 

TM !56:"li.;2.M* *;577,.13,-.12e5,-.°03067, 

j=z 

IF( Z . LE . 37 . ) THEM 
CCOEF=C( 1 ,<J) 

RETURN 

RETURN 

ELSE 

DO 3 1 1 ■= 1 , 6 

IIZ=II . 

IF(J.LT.ZB( 11) ) GO TO 4 

3 CONTINUE 

4 CCOEF=CC (I,1IZ-1)+ (CC ( I , I IZ ) -CC (I,IIZ-1))* 

-(Z-ZB(IIZ-l) ) /( ZB( I IZ ) -ZB( 1 IZ- 1 ) ) 

RETURN 
END IF 
END 

****************************************************** 

SUBROUT I NE MGAUSS ( A , B , N , SUM , FUNC , FOFX , NUMBER ) 

cM.nuniiTtNF MGAUSS USES A GAUSS- LEGENDRE QUADRATURE FORMULA TO 
PERFORM A 5 POINTS PER SUBINTERVAL NUMERICAL INTEGRATION 

A= LOWER LIMIT 
B=UPPER LIMIT 
N=NUMBER OF INTEftVALS 
«?!fM=RESULT Of INTEGRATION 

fUNONAME Of SUBROUTINE THAT GIVES INTEGRAND 

FOFX* INTEGRAND ^ _ lt , 

NUMBER=NUMBER OF INTEGRALS BEING DONE 

DIMENSION U( 5 ) , R ( 5 ) ,SUM( 1) ,FOFX( 1) 

DO 10 LL=1, NUMBER 
10 SUM(LL)=0 .0 

IF(A.EQ.B) RETURN 

u(i) Are the zeros of legendre polynomial of deg. 5 


U( 1)=. 42556283050918 



nooanfj nno 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 


c 

c 


fc.2fi330230293537 
s.160295215850480 
a. 067468316655508 
-.013046735741414 

R(l) ARB THE WEIGHTS SUCH THAT THE 



INTEGRATION IS EXACT 


Rl 1)-. 147762112357376 
r(2)-. 13463335965499 
R(3)»- 109543181257991 
R( 4 )“• 07 472 56 74575290 
R(5)=. 033335672154344 
FINE-N 

DELTA=FlNR/(fi-A) 


FIND INITIAL VALUE FOR INTERVAL 


DO 30 K=1,N 
XI =K- 1 

FiNE-A+XI/DELTA 


EVALUATE FUNCTION AT 5 POINTS PER INTERVAL 


DO 20 11-1,5 

UU-U(II)/DELTA+FINE 
CALL FUNC ( UU , FOFX ) 

ADD TO SUM FOR EACH FUNCTION 

DO 20 NN=1, NUMBER 
20 SUM(NN)=R( II ) *FOFX(NN)+SUM(NN) 

DO 30 11-1,5 

UU=(1.0-U(II))/DELTA+FINE 

CALL FUNC (UU,FOFX) 


DO 30 NN=1, NUMBER 
30 SUM ( NN ) =R ( 1 1 ) * FOFX ( NN ) +SUM ( NN ) 


DO 40 1=1, NUMBER 
40 SUM( I ) =SUM( I ) /DELTA 
RETURN 
END 

* * * * SUBROUT I NE IUNI ( NMAX , N , X 


NTAB, Y , IORDBR ,X0,¥0 , IPT, I ERR) ******** 


SUBROUT I NE I UN I ( NM AX , N , X , NT AB , Y , I ORDER , XO , YO , I PT , I ERR ) 

SUBROUTINE IUNI USES A 1-ST OR 2-ND ORDER LAGRAHG1AN INTERPOLATION 
FOR A MONOVARIATE FUNCTION Y=F(X) 


TAKEN FROM LANGLEY MATH. LIBRARY 

DIMENSION X(NMAX) ,Y( NMAX, NTAB) , YO(WTAB) 

NM1-N-1 
IERR-0 

J=i 

DELX=X ( 2 ) -X ( 1 ) 

IF ( IORDER • EQ. 0) GO TO 10 
IF (N.LT. 2) GO TO 20 
GO TO 50 
10 IERR--1 
GO TO 30 
20 IERR--2 
C 



non 


Jo DO 40 NT“1 »NtAB 
¥0(M*)*¥( l*Nt ) 

4d coNtINUft 

c 

RfeTtlRM 

50 IF (IPT .GT. -1) GO TO 65 
IF (DELX . EQ. 0) GO TO 190 
IF (N .EQ. 2) GO TO 65 
C 

DO 60 J=2,NMl 

IF (DELX * (X(J+1 )-X( J) ) ) 190,190,60 
60 CONTINUE 
C 

65 IF (IPT .LT. 1) 1PT=1 

IF (IPT -GT. NMl ) IPT=NM1 
IN= SIGN ( 1 -O f DELX *( XO-X(IPT))) 

70 P= X(IPT) - XO 

IF (P* ( X ( IPT +1)- XO)) 90,180,80 
80 IPT =IPT + IN 

IF (IPT.GT.O -AND. IPT -LT. N) GO TO 70 
IERR=-4 
IPT=IPT- IN 

90 IF (IORDER -GT. 1) GO TO 120 
DO 100 NT“l f NTAB 

YO ( NT ) =¥ ( I PT , NT ) + ( ( Y ( IPT+ 1 , NT) - Y ( IPT , NT) ) * ( X0-X( IPT) ) )/ 

-(X( IPT+1 ) -X( IPT) ) 

100 CONTINUE 

IF ( 1ERR -EQ. -4) IPT=IPT+IN 
RETURN 

120 IF (N .EQ. 2) GO TO 200 

IF (IPT .EQ. NMl) GO TO 140 
IF (IPT .EQ. 1) GO TO 130 

IF (DELX * ( XO-X ( IPT- 1 ) ) . LT. DELX* ( X( IPT+2 ) -XO ) ) GOTO 140 
130 L“IPT 

GO TO 150 
140 L“IPT -1 
150 Vl=X(L) -XO 
V2=X(L+1 ) -XO 
V3“X ( L+2 ) -XO 

DO 160 NT“1,NTAB 

¥Yl=( Y ( L, NT ) * V2 - Y( L+l , NT) * Vl)/(X(L+1) - X(L)) 
YY2“(Y(L+1,NT)*V3-Y(L+2,NT) *V2 ) /(X ( L+2 ) -X( L+l ) ) 

160 Y0(NT)=( YYl*V3 -YY2*Vl)/(X(L+2)-X(L) ) 

IF ( IERR .EQ. -4) IPT=IPT + IN 
RETURN 

180 IF( P .NE. 0) IPT=IPT +1 

DO 185 NT=1,NTAB 
Y0(NT)“Y( IPT, NT) 

185 CONTINUE 

RETURN 

190 IERR=J +1 
RETURN 
200 IERR=-3 
RETURN 
END 

***** SUBROUT I NE IBI ( NX , X , NY , Y , MAXF, F, IORDER, IPT,X0,Y0,Z, IERR ) ****** 
SUBROUT INE IBI(NX,X,NY,Y, MAX F , F , IORDER , I PT , XO , YO , 2 , I ERR ) 



nnnn n n a 


c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


c 


SUBROUTINE 181 USES A 1-ST OR 2-ND ORDER LAGRANGIAN INTERPOLATION 
frott A BlVAftlAtfi PUNcttoN 2*P(X,Y) 


tAkEN PRoM LANGLEY MAtH LIBRARY 

REQUIRED ROUTINES QXZ 114 , QXZ 115, QXZ 116, QXZ 117 

DIMENSION X(NX) ,Y(NY) f P(MAXP,NY) , P( 3 ) , S( 3 ) , 1PT( 2 ) , IORDER( 2 ) 

I F ( I PT ( 1 ) . NE . - 1 ) GO TO 10 
IERR = - 1 

IF(QXZ 114 (NX , X( 1) ,IPT(1) ) .LT.O) RETURN 
IF(QXZll4(NY,Y< 1) , IPT(2) ) .LT.O) RETURN 
IPT ( 1 ) = 1 
IPT( 2 ) = 1 
IPT ( 1 ) = IPT ( 2 ) = 1 
10 IERR s -4 

I F I IORDER ( 1 1 . NE . 1 . AND . IORDER ( 1 ) . NB . 2 ) RETURN 
IF( IORDER( 2 ) . NE . 1 • AND . IORDER ( 2 ) . NE . 2 ) RETURN 
lERR 9 -2 

IFINX-l . LT. IORDER ( 1 ) ) RETURN 
IF(NY-1 . LT. IORDER (2 ) ) RETURN 
IERR-0 

CALL QXZ 116(X0,X,IPT( 1) , NX ,M , I ERR ) 

CALL QXZll6(Y0,Y,IPT(2) , NY , N, I ERR ) 

CALL QXZll5(IPT(2), IORDER ( 2 ) , Y0 , KL, KR , Y , NY ) 

RM=KL+ 1 

CALL QXZll7(YO,Y(KR),Y(KM),Y(KL),P,IORDER(2)) 

CALL QXZ 1 15(1 PT(1), IORDER ( 1 ) » X0 , LL, LR , X, NX ) 

LM=LL+1 

CALL QXZ 1 17 ( XO , X ( LR) , X ( LM ) , X ( LL) , S , IORDER ( 1 ) ) 

Z=0. 

LP1=0 


DO 30 I=LL, LR 
LPl=LPl+l 
SUM=0 . 

KPl=0 

DO 20 J=RL,kR 
KPl=KP 1+ 1 

20 SUM=P(KP1 ) *F( I , J)+SUM 


30 Z=SUM*S ( LP 1 ) +Z 

IF (IBRR.GE.O) RETURN 

IPT( 1 ) =M 

IPT(2)=N 

RETURN 

END 


************** **FUNCTION QXZl 14 (M,VAR, IPOS) ********************** 
FUNCTION QXZll4(M, VAR, IPOS) 

FUNCTION QXZl 14 CHECKS TO SEE IF SEQUENCE VAR IS INCREASING 
TAKEN FROM LANGLEY MATH. LIBRARY 

DIMENSION VAR(M) 

QXZl 14=1 
K=M- 1 

DO 10 L=1 , K 
N=L+ 1 

IF(VAR(N|-VAR(L) .GT.O. )GO TO 10 

QXZ114=-1 

IPOS=N 



nnooon nnn 


io coMttMUB 

ftfetURM 
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C 

C 

c 

c 

c 

c 


C**4********SUBROUTINE QXZ 1 15 (IPT,IORDER,S,IL,IR,VAR,N) ************* 
C SUBROUTINE QX2l i5( IPT, loRDER, S, IL, IR, VAR,N) 

SUBROUTINE QX2115 DfeCIDES THE INDICIES OF X,Y,2 ARRAY PASSES 
TO ARRAY VAR 

TAKEN PROM LANGLEY MATH. LIBRARY 

DIMENSION VAR(N) 

1 F ( S • LE . VAR (2)1 GO TO 10 
IF(S.GE.VAR(N-1) )GO TO 20 
IL*lPT 

IFCS.LE.VAR(IPT) ) IL=IL-1 

SlSSS-vSSi u? *«* ' ' ' **-«- 1 

GO TO 30 
10 IL*1 

GO TO 30 
20 IL=N-IORDER 
30 IR=IL+IORDER 
RETURN 
END 

*********** SUBROUT I NE QXZ 1 1 6 ( S , VAR , I POS , NBND , I F , I ERR )************ 


SUBROUTINE QXZ 1 1 6 ( S , VAR , I POS , NBND , I F , I ERR ) 

SUBROUTINE QXZll 6 DECIDES THE INDEX OF X,Y,Z ARRAY ( VAR ) 
NEAREST TO X 0 ,Y 0 ,Z 0 

taken From langley math, library 


DIMENSION VAR (NBND) 
LOGICAL ABOVE, BELOW 


C 


C 


ABOVE*. FALSE. 

BELOW* . FALSE . 
ABOVE=BELOW* .FALSE . 

I F ( I POS . EQ . 0 ) 1 POS * 1 
10 IF(S-VAR( IPOS) ) 20 , 50 , 
20 ABOVE*. TRUE. 


30 


IUP=IPOS 

IF( ABOVE. AND. BE LOW) GO TO 

IPOS=IPOS-l 

I F ( I POS . GE . 1 ) GO TO 10 

IERR =-3 

IPOS=l 


IF =0 

RETURN 

30 BELOW*. TRUE. 

LOW* I POS 

IF (ABOVE. AND. BE LOW) GO TO 
IPOS*IPOS+l 

IF( IPOS • LE . NBND) GO TO 10 
IERR =-3 

if*nbnd 

IPOS*NBND 

RETURN 


40 


40 



40 lPOS=LOW 

IF(ABS(S-VAft( 1UI 5 ) ) .Lf.ABS(S-VAR(LoW) ) )tt>OS=lUP 
50 1F-IPOS 
RETURN 
END 
C 

C* ********* * SUBROUt INE QXZll7(S,VlRl, VAR 2 , VAR3 , Q , IORDBR )************ 

C 

SUBROUTINE QXZ 1 1 7 ( S , VAR 1 , VAR2 f VAR3 ,Q, IORDER ) 

C 

C SUBROUTINE QXZ 113 CALCULATES THE LAGRANGIAN COEFF. FOR 
C 1-ST OR 2-ND ORDER INTERPOLATION 
C 

c TAkEN FROM LANGLEY MATH. LIBRARY 

C 

DIMENSION Q( 3 ) 

C 

IF( IORDER. EQ. 2 )GO TO 10 
Tl=VAR3-VARl 
T2=S-VAR1 
T3=S-VAR3 
Q( 1 )=T2/Tl 
Q(2)=-T3/Tl 
RETURN 
10 T1-S-VAR2 
T2=S-VAR1 
T3=*S-VAR3 
Q( 1 )=Tl*T2 
Q(2)=T3*T2 
Q( 3 ) =T3*T 1 
T1=*VAR3-VAR2 
T2-VAR3-VAR1 
T3=VAR2-VARl 
Q(1)=Q(1)/(T1*T2) 

Q( 2 ) =Q( 2 )/(-Tl*T3) 

Q( 3 ) =Q( 3)/(T2*T3) 

RETURN 

END 

C 

Q* ************************************************ 

C 

FUNCTION TEXP(X) 

C 

IF(X.LT.-88. ) THEN 
TEXP=EXP ( -88 . ) 

RETURN 

ELSE IF(X.GT. 88 . ) THEN 
TEXP=EXP( 88 . ) 

RETURN 

ELSE 

TEXP=EXP ( X ) 

END IF 
RETURN 
END 
C 

c 

SUBROUTINE SUTPAR ( NX , NY , X , Y , NZ , Z , IBNDSW , BNDX 1 , ENDXN , ENDY 1 , 

$ ENDYN , ENDXY, SIGMA , M , XX , YY , IW, ZZ , WK , 

$ IERR ) 

C 

C SUBROUTINE SUTPAR INTERPOLATE A BIVARIATE SPLINE UNDER TENSION 
C AT A VECTOR OF POINTS, RETURNING THE FUNCTION VALUE, THE TWO FIRST 
C PARTIAL DERIVATIVES, AND THE THREE SECOND PARTIAL DERIVATIVES 

C 

C TAKEN FROM LANGLEY MATH. LIBRARY 



c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


ttfeQUittfeb RoOtlMBS - QXZ060,QXZ063,QX2062,QXZll9,QXZlll,QX206l,QXZll2, 

FORMAL PARAMEtERS 

INtEGER t ENDS# (8) , 1 ERR , I W , M , NX , NY , NZ 

BEAL ENDX 1 ( MY ) , ENDXN( NY ) , ENDY 1 (NX ) , ENDYN(NX ) , ENDXY ( 4 ) , SIGMA, WK( * ) 
REAL X(NX) ,XX(M),Y(NY),YY(M) ,Z(NZ,NY) ,ZZ(6,M) 

INTERNAL VARIABLES 

INTEGER INDEX ,K 

REAL ABSIGMA, DELX , DELY , ONE , SIGMAX , SIGMAY , ZERO 
INITIALIZE CONSTANTS 

DATA ZERO, ONE/O. 0E0, 1 .0E0/ 

INITIALIZE IERR AND CHECK ERROR RETURNS 


IERR =* 1 

IF (NX .LE. 1) GO TO 9000 
IF(M.LE.O) GO TO 9000 
IF (NY .LE. 1) GO TO 9000 
IERR s 2 

DELX - ( X ( NX ) - X ( 1 ) ) /( NX - 1) 
IF (DELX .LE. ZERO) GO TO 9000 
DELY - ( Y ( NY ) - Y(1))/(NY - 1) 
IF (DELY .LE. ZERO) GO TO 9000 
IERR = 0 


DENORMALIZE THE TENSION FACTOR SIGMA FOR X AND * 
IF SIGMA IS LESS THAN l.E-30, IT IS REDEFINED TO ZERO 
(WITH NO SIGNIFICANT EFFECT) TO AVOID UNDERFLOW. 


ABSIGMA = ABS( SIGMA) 

IF (ABSIGMA .LT. l.E-30) 
SIGMAX = ABSIGMA / DELX 
SIGMAY = ABSIGMA / DELY 


ABSIGMA =0.0 


BRANCH IP QXZ 119 HAS ALf 

IF (1W .EQ. 1) GO TO 30 
INDEX = NX*NY*3 + 1 
IF ( IENDSW( 3 ) .NE. 2) G( 
IENDSW(5) = 0 
ENDXY(l) = ZERO 
IENDSW( 6 ) = 0 
ENDXY ( 2 ) = ZERO 
10 IF (IENDSW(4) .NE. 2) G< 
1BNDSW( 7 ) = 0 
ENDXY ( 3 ) = ZERO 


TO 10 


ENDXY (1),ENDXY(ZJ, 
DELY , SIGMAY r 0 , 0 , WP 
t (IERR .NE. 0) GO TO 9000 

MUST BE SET tO 1 FOR OUTPUT 


C 


IW = 1 
30 CONTINUE 



non nnn nnn nnn non non non n nnnnnnnnn non non 


BEGIN Loot* Fott IMfERFoLATlON AMD DERIVATIVE EVALUATION 
DO 90 R=1,M 

CALL QXZl ll(XX(k) , YY(k) ,ZZ( 1 ,k) ,ZZ(2,K) ,ZZ(3,K) ,ZZ(5,k) , 

* ZZ(4,k),ZZ(6,K),NX,NY,X,Y, Z ,NZ ,WK, SIGMA, I ERR ) 

90 CONTINUE 
RETURN 

00 IF(IERR .EQ. 1) IERR - -2 
IF( IERR .EQ. 2) IERR = -1 
RETURN 
END 




SUBROUTINE QXZ060 ( DELI , DEL2 , SIGMA, C 1 ,C2 , C3 , N) 

SUBROUT1ME QXZ060 DETERMIME COEFFICIENTS Cl , C2 , AND C3 
USED TO DETERMIME ENDPOINT SLOPES. GIVEN THE 
POINTS ( Xl , Yl ) , ( X2 , Y2 ) , AND (X3,Y3) , 

THE ENDPOINT SLOPE IS Cl*Yl + C2*Y2 4- C3*Y3 * 

IF N > 2 AND Cl *Y 1 + C2*Y2 IF N =2. 

FORMAL PARAMETERS 

INTEGER N 

REAL Cl, C2,C3, DELl, DEL2, SIGMA 
INTERNAL VARIABLES 

REAL COSHM1 , DEL, DENOM, ONE, SIN, TWO, ZERO 
INITIALIZE CONSTANTS 

DATA ZERO, ONE, TWO/O.OEO, 1 .0E0,2.0E0/ 

TEST WHETHER N ■> 2 
IF (N .EQ. 2) GO TO 20 

IF (ABS(SIGMA*DEL2) .GT. 740.0) GO TO 30 

INITIALIZE THE VARIABLE DEL 

DEL = DEL2 - DELl 

TEST WHETHER SIGMA = ZERO 

IF (SIGMA .ME. ZERO) GO TO 10 

THIS CODE IS FOR A STANDARD CUBIC SPLINE 

Cl *> -(DELl+DEL2)/(DELl*DEL2) 

C2 = DEL2/( DELl *DEL) 

C3 = -DELl/(DEL2*DEL) 

GO TO 9000 

THIS CODE IS For A SPLINE UNDER TENSION 

10 CALL QXZ062(DUMMY,COSHM1,S1GMA*DBL1,1) 

CALL QXZ062 ( DUMMY, COSHM2 , SIGMA*DEL2 , 1 ) 

SIN - TWO * SINH ( SIGMA* ( DEL2+DBL1 ) /TWO) * SINH( SIGMA*DBL/TWO) 
DENOM - ONE / (COSHMl*DEL - DELl*SIN) 

Cl => SIN * DENOM 
C2 = -COSHM2 * DENOM 
C3 « COSHM1 * DENOM 



GO TO 9000 
C 

c to Avoid overflow if S1gMA*deL 2 IS Too LARGE 
c 

30 c3 • o.o 
c 

c TWO COEFFICIENTS 

c 

20 Cl - -ONE/DELl 
C2 = -Cl 
9000 RETURN 
END 
C 

Q *****A*********************************************** 

C 

SUBROUT IRE QXZ06 1 ( DI AG, SDIAG, SIGMA, DEL) 

C 

C SUBROUTINE QXZ061 COMPUTE THE DIAGONAL AND SUPERDIAGONAL TERMS 
C OF THE TRIDIAGONAL LINEAR SYSTEM ASSOCIATED 
C WITH SPLINE UNDER TENSION INTERPOLATION. 

C 

C FORMAL PARAMETERS 

C 

REAL DEL,DlAG, SDIAG, SIGMA 
C 

C INTERNAL VARIABLES 
C 

REAL COSHM, DENoM, SIGDEL, SINHM 
REAL* 8 SIXTH , THIRD 
C 

C INITIALIZE CONSTANTS 
C 

DATA ZERO, SIXTH, THIRD/ 

$0.0E0, ' 17155252525252525253 '0, ' 17165252525252525253*0/ 

C 

C TEST WHETHER SlGMA IS ZERO 
C 

IF (SIGMA .NE. ZERO) GO TO 10 
C 

C THIS CODE IS FOR THE STANDARD CUBIC SPLINE 
C 

DlAG - DEL*THIRD 
SDIAG - DEL*SIXTH 
GO TO 9000 
C 

C THIS CODE ig FOR A SPLINE UNDER NON-ZERO TENSION 
C 

10 SIGDEL ** SIGMA*DEL 

CALL QXZ062 (SINHM, COSHM, SIGDEL, 0) 

DENOM - SIGMA* SIGDEL* ( 1 . +SINHM ) 

DIAG ■ (COSHM- SINHM) /DENOM 
SDIAG ” SINHM/DENOM 
C 

9000 RETURN 
END 
C 

Q *************** *****4 ******************** *********** 

C 

SUBROUTINE QXZ062 ( SINHM, COSHM, X, ISW) 

C 

C SUBROUTINE QXZ062 RETURNS APPROXIMATIONS TO 
C SINHM(X) =SINH( X ) /X- 1 

C COSHM ( X ) "COSH ( X ) - 1 

C COSHMM(X)=(COSH(X)-l-X*X/2)/(X*X) 

C WITHOUT RELATIVE ERROR LESS THAN 4.0B-14. 

C 



no non non 


INTEGER isw 

REAL SINHM,COSttM,X 

DATA SP14/.227581660976348E-7/, 

* SP13/.612189863171694E-5/, 

* SP12/.715314759211209E-3/, 

* SP11/.398088289992973E-1/, 

* SQ12/. 2063827014 13725E-3/ , 

* SQl 1/- . 6 1 14 70260009508E- 1 / , 

* SQlO/. 599999999999986E+1/ 

DATA SP25/. 129094 15B037272E-9/ , 

* SP24/.473731823101666E-7/, 

* SP23/.849213455598455E-5/, 

* SP 22/. 833264803327242 E-3/ f 

* SP21/.425024142813226E-1 / , 

* SQ22/. 106008515744821E-3/, 

* SQ21/-.449855169512505E-1/ , 

* SQ20/ . 6000000002686 19E+ 1/ 

DATA SP35/. 155193945864942B-9/, 

* SP34/.511529451668737E-7/, 

* SP33/.884775635776784E-5/, 

* SP32/.850447617691392B-3/, 

* SP31/.428888148791777B-1/, 

* SQ32/.933128831061610E-4/, 

* SQ3 1/- .4 2667 7 570538 50 7E-1/ , 

* SQ30/.600000145086489E+1/ 

DATA SP45/. 188070632058331 E- 9/ , 

* SP44/.545792817714192E-7/, 

* SP43/.920119535795222E-5/ , 

* SP 42/. 86655939 167 2985E-3/ , 

* SP41/.432535234960858B-1/, 

* SQ42/.824891748820670E-4/, 

* 5Q4 1/-. 404938841 67 2262E-1/ , 

* SQ40/. 600005006 2 838 34 E+ 1/ 

DAtA CP5/.552200614584744B-9/, 

* CP4/.181666923620944E-6/, 

* CP3/.270540125846525B-4/, 

* CP2/ . 2062707 19503934E-2/ , 

* CP 1/. 74443720556904 0E-l/ f 

* CQ2/.514609638642689B-4/, 

* CQl/-. 177792255528382E-1/, 

* CQO/ . 200000000000000E4 1/ 

DATA ZP4/.664418805876835E-8/, 

* ZP3/.218274535686385E-5/, 

* ZP2/.324851059327161E-3/, 

* ZP1/.244515150174258E-1/, 

* ZQ2/. 6 16 16578230662 IE- 3/ , 

* ZQl/- . 2 1 3 1636395794 25E0/ , 

* ZQ0/.240000000000000E+2/ 

AX - ABS(X) 

IF ( ISW .GE. 0) GO TO 5 

SINUM APPROXIMATION 

IF (AX .GT. 3.9) GO TO 2 
XS - AX* AX 

IF (AX .GT. 2.2) GO TO 1 

SINUM APPROXIMATION ON (0.,2.2) 

SINHM - XS*( ( ( (SPl4*XS+SPl3)*XS+SPl2) *XS+SPl 1)*XS+1. )/ 
( ( S9l2*XS+SQl 1 ) *XS+SQ10) 

RETURN 

SINUM APPROXIMATION ON (2. 2 ,3. 9) 



oooon n n o non non non non 


c 


1 SINMM * XS*( ( ( ((SP25*XS4SP24)*XS4SP23)*XS4SP22)*XS4SP2l) 

*XS41. )/( (SQ22*XS4SQ2l)*XS4SQ20) 

RETURN 

2 IF (AX .Gt. 5.1) GO TO 3 
SINMM APPROXIMATION ON (3. 9, 5.1) 

XS = AX*AX 

SIMHM = XS*( ( ( ( (SP35+XS4SP34)+XS4SP33)*XS4SP32)*XS4SP31) 
♦XS+1 . )/( (SQ32*XS4SQ31 )*XS4SQ30) 

RETURN 

3 IF (AX .GT. 6.1) GO TO 4 
siNMM Approximation on (5.i,6.i) 

xs - Ax*Ax 

SINMM - XS*( ( ( ( (SP45*XS+SP44 ) *XS+SP43)*XS+SP42)*XS+SP41) 
*XS+1. )/( ( SQ4 2*XS+SQ4 1 ) *XS4SQ40 ) 

return 

SINMM APPROXIMATION ABOVE 6.1 

4 EXPX = EXP(AX) 

SINMM » (EXPX-1./RXPX)/(AX4AX)-1. 

RETURN 

COSMM AND (POSSIBLY) SINMM APPROXIMATION 

5 IF (ISW .GB. 2) GO TO 7 
IF (AX .GT. 2.2) GO TO 6 
XS ** AX*AX 

COSMM ** XS*( ( ( ( (CP5*XS4CP4 ) *XS+CP3 ) *XS+CP2 ) *XS+CPl ) 

*XS4l . )/( (CQ2*XS4CQl )*XS4CQ0) 

IF (ISW .EQ. 0) SINMM « XS* ( ( ( ( SP14 *XS4SP1 3 ) *XS4SPl2 ) 
♦XS4SP11 ) *XS4l . )/( (SQl2*XS4SQll)*XS4SQlO) 

RETURN 

6 EXPX » EXP(AX) 

COSMM = ( EXPX+ 1 . /EXPX ) /2 . - 1 . 

IF (ISW .EQ. 0) SINMM - ( EXPX- 1 . /EXPX ) /(AX4AX )- 1 . 

RETURN 

COSMMM AND (POSSIBLY) SINMM APPROXIMATION 

7 XS = AX* AX 

IF (AX .GT. 2.2) GO TO 8 

COSMM «* XS*( ( ( (ZP4*XS4ZP3)*XS4ZP2)*XS4ZPi)*XS41. )/ 

( (ZQ2*XS4ZQl ) *XS42Q0) 

IF (ISW -EQ. 3) SINMM - XS* ( ( ( ( SP14 *XS4SPl 3 ) *XS4SPl2 ) 
*XS4SP 1 1 ) *XS+ 1 . )/( (SQ12*XS4SQ1I)*XS4SQ10) 

RETURN 

8 EXPX - EXP ( AX ) 

COSMM - ( ( EXPX4 1 . /EXPX-XS J/2.-1 . )/XS 

IF (ISW .EQ. 3) SINMM = (EXPX-1 ./EXPX)/(AX4AX)-1 . 

RETURN 

END 

**A+#***********4******************#*#4***«***** 


INTEGER FUNCTION QXZ063 (T,X f N) 

FUNCTION QXZ063 DETERMINES THE INDEX OF THE 1NT 
(DETERMINED) BY A GIVEN INCREASING SEQUENCE) 
WHICH A GIVEN VALUE LIES. 


INTEGER N 



nnnonn 


REAL T,X(N) 

C 

save l 
data i /l/ 
c 

TT = T 

CHECK FOK ILLEGAL 1 

It - (I .GE. N) I - N/2 

CHECK OLD INTERVAL AND EXTREMES 

It- (TT .LT. X( I ) ) THEN 
It- (TT .LB. X ( 2 ) ) THEN 
I - 1 

QXZ063 <= 1 
RETURN 
ELSE 
IL = 2 
111 - I 
END IF 

ELSE IF (TT .LB. X(l+1)) THEN 
QXZ063 - I 
RETURN 

ELSE IF (TT .GB. X(N-1 ) ) THEN 
I = N- 1 
QXZ063 = N- 1 
RETURN 
ELSE 
IL - 1+1 
111 - N-l 
END IF 

BINARY SEARCH LOOP 

1 I - (IL+IHJ/2 

IF (TT .LT. X ( I ) ) THEN 
IH = I 

ELSE IF (TT .GT. X(I+ll) THEN 
IL - 1+1 
ELSE 

QXZ063 « I 
RETURN 
END IF 
GO TO 1 
END 


************* *************************************** 

SUBROUTINE QXZlll (XX,Y¥,ZZ,ZX,Z¥, ZXX , ZX¥,Z¥Y,M,N,X,Y, 
* Z, IZ,ZP, SIGMA, IBRR) 

THIS SUBROUTINE EVALUATES THE FUNCTION VALUE. THE TWO 

DERIVATIVES, AND THE THREE SECOND PARTIAL 
DERIVATIVES OF A TENSOR PRODUCT SPLINE UNDER TENSION IN 

^t^B 1ABLeS * ™ E SUBROUTINE QXZ119 SHOULD BE CALLED 
EARLIER TO DETERMINE CERTAIN NECESSARY PARAMETERS. 

INTEGER M,N, IZ 
INTEGER QXZ063 

RBM, CC(10MJ(2),XJ(2),XXJ(2|,YYJ(2),XYYJ(2),XXYYJ(2| 




DE NORMALIZE tENSloN FACTOR IN X AND Y DIRECTIONS 

ABsIgMA - AflS( SIGMA) 

IF(ABSIGMA .LT. l.B-30) ABSIGMA => 0.0 
SIGMAX - ABSIGMA*PLOAT(M-l)/(X(M)-X(l)) 

SIGMAY - ABSIGMA*FLOAT(N-l)/(Y(N)-Y( 1) ) 

DETERMINE X INTERVAL 

IM1 - QXZ063(XX,X( 1) ,M) 

I = IM1 + 1 

determine y interval 

JM1 - QXZ063(YY,Y( 1) ,N) 

J - JM1+1 

OBTAIN X MESH DISTANCES 

DELI - XX-X(IMl) 

DEL2 - X(I)-XX 
DELS - X ( I ) -X ( IM 1 ) 

SET ERROR CODE TO PREVENT OVERFLOW DUE TO EXCESSIVE EXTRAPOLATION. 

repp s 4 

IF C SIGMAX*ABS( DELI ) .GT. 650.0 .OR. SIGMAX* ABS( DBL2 ) .GT. 650.0) 
* ' GO TO 9000 

I ERR - 0 

OBTAIN COMMON COEFFICIENTS 


CC(1) “ DEL2/DELS 
CC( 2 ) = DELl/DELS 
D - - 1 . /DELS 
CC( 5 ) =* D 
CC( 6) - -D 

IF (SIGMAX .NE. 0.) GO TO 1 
OBTAIN COEFFICIENTS FOR CUBIC SPLINE 


CON = -DELI *DEL2/( 6 . *DELS ) 

CC( 3 ) - ( DEL2+DELS ) *CON 
CC( 4 ) = ( DELl+DELS ) *CON 

CC (7j = - ( 2 . *DEL2 *DEL2-DELl * (DEL2+DELS ) ) /( 6 . *DELS) 
CC(8) - ( 2 . *DELl *DELl-DEL2* ( DELl+DELS) ) /( 6 . *DELS) 
CC ( 9 ) - CC(1) 

CC(10) = CC ( 2 ) 

GO TO 2 

OBTAIN COEFFICIENTS FOR SPLINE UNDER TENSION 


1 SIGDEL » SIGMAX* DELS 

CALL OXZ062 ( SINHMl ,COSHMl , SIGMAX*DELl , 0) 
CALL QXZ062 ( SINHM2 ,COSHM2 , SIGMAX*DEL2 , 0) 
CALL QXZ062(SINHMS, DUMMY , SIGDEL ,-l) 

SINHS - 1. + SINHMS 
CON - 1 ./( SIGMAX* SIGDEL* SINHS) 

CC( 3 ) = CON*DEL2 * ( S1NHM2 -SINHMS ) 

CC(4) - CON*DELl * ( SINHMl -SINHMS) 

CC(7) - -CON* (COSHM2-SINHMS) 

CC( 8 ) - CON* (COSHM 1- SINHMS ) 

CC(9) - ( DEL2* ( 1 • +SINHM2 ) ) /( DELS* SINHS ) 
CC(IO) - ( DELl * ( 1 . +SINHM I ) ) /( DELS* SINHS) 

INTERPOLATE IN X DIRECTION 



non nnn nno 


i CALL QX4112 
I t i ) } 

CALL QXil 12 (CC*ZJlMl,J) f ZP(tMl^,2),ZJ(2),XJ(2) f 

*CALL QXZ112 (CcJit>(lMl,JMl, 1) , ZP( 1M1 , JM1 , 3 ) ,YYJ( 1) , 

* XYY J ( 1) , XXYYJ ( 1) ) 

CALL QXZ112 (CC,ZP(lMl,J,l),ZP(lMl,J,3),YYJ(2), 

* XYYJ (2) ,XXYYJ(2) ) 

OBTAIN Y MEStl DISTANCES 

DELI - YY-Y(JMl) 

DEL2 - Y(J)-YY 
DELS - Y(J)-Y(JM1) 

SET ERROR CODE TO PREVENT OVERFLOW DUE TO EXCESSIVE EXTRAPOLATION. 

TF R (SIGMAY*ABS(DELl) .GT. 650.0 .OR. SIGMAY*ABS( DEL2 ) .GT. 650.0) 
* ' GO TO 9000 

IERR = 0 

OBTAIN COMMON COEFFICIENTS 

CC(1) - DEL2/DELS 
CC ( 2 ) - DELl/DELS 
D » - 1 . /DELS 
CC( 5 ) - D 
CC( 6 ) = -D 

IF (SIGMAY .ME. 0.) GO TO 3 
OBTAIN COEFFICIENTS FOR CUBIC SPLINE 


CON - 
CC( 3 ) 
CC(4) 
CC( 7 ) 
CC(8) 
CC( 9 ) 


-DELl *DEL2/( 6 . *DELS) 

= (DEL2+DELS) *CON 

- (DELl+DELS) *CON 

= - ( 2 . *DEL2 *DEL2-DELl * (DEL2+DELS ) )/(6.*DBLS 
. (2.*DELl*DELl-DEL2*(DELl+DELS) )/(6.*DELS) 

- CC(1) 


CC (10) - CC ( 2 ) 

GO TO 4 

OBTAIN COEFFICIENTS FOR SPLINE UNDER TENSION 

3 SIGDEL « SIGMAY *DELS 

CALL OXZ062 ( SINHM1 ,C0SI1M1 , SIGMAY*DBLl , < 

CALL QXZ062 ( SINHM2 ,COSHM2 , SIGMAY*DEL2 , ( 

CALL QXZ062(SINHMS, DUMMY , SIGDEL 

SINUS - 1 . +SINHMS 

CON = 1 ./( SIGMAY* SIGDEL* SINHS ) 

CC(3) = CON*DEL2* ( SINHM2-SINHMS) 

CC(4) = CON*DELl* (SINHM1-SINHMS) 

C C( 7 ) - -CON* (COSHM2-S1NHMS ) 

CC(8 ) » CON* (COSHM1-SINHMS) 

CC(9) - ( DEL2* ( 1 . +SINHM2 ) ) /(DELS* SINHS) 
CC(10) - ( DELl * ( 1 . +SINHM1 ) )/(DELS*SlNHS) 


’ INTERPOLATE IN Y DIRECTION 
1 

4 CALL QXZ112 (CC ,ZJ,YY<J,ZZ f ZY r ZYY) 

CALL QXZ112 (CC, XJ , XYYJ, ZX, ZXY f DUMMY ) 

ZXX = CC( 1)*XXJ( 1)+CC(2)*XXJ(2)+CC(3)*XXYYJ(1) + 
* CC( 4 ) *XXYYJ ( 2 ) 

9000 RETURN 
END 



non 


c 

c 

c 

c 

c 


SUBROUTINE QXZll2 (A,Y,YP,YY,YYP,YYPP) 
REAL A(l0),Y(2),YP(2) 

YY - A( 1 ) *Y( 1 


1 . A ^ . a . . 4 > 4 n k . trtk < 1 I 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


yyp «s A(S 

YYPP - A(9)*YP(1)+A(10)*YP(2) 

return 

END 

♦ 4444********* A***************************************”*** 

«M»on M w»t igfg*; • 

| ZP, TEMP, IERR) 

SUBROUTINE QXZ119 COMPUTE PARAMETERS NECESSARY TO BE ABLE TO 
INTERPOLATE POINTS ON A BIVARIATE SPLINE 
FUNCTION UNDER TENSION PASSING THROUGH A 
RECTANGULAR GRID OF FUNCTIONAL VALUES. SEE 
STIBI (El.X). 

formal parameters 

INTEGER IERR, IOPTX, IOPTY ,ISLPSW(8),IZ,M,N 

REAL DELX DELY,SIGMAX,SIGMAY,TEMP(*) ,X(M) ,Y(N) ,Z(IZ,N) ,ZP(M,N,3) 

SSE ZXH^ N) ,Zxi (H) , ZYH(H) , ZY1 (M| , ZXYMN,ZXYH1 *ZXY IH, ZXY1 1 


c 

c 

c 

c 

c 

c 


c 

c 

c 


INTERNAL PARAMETERS 

INTEGER I , I BAR , 1BAKP 1 , IMl , IPl ,J , JBAK , JBAKPl , JM1 , JPl ,MM1 ,MPl ,NM1 
INTEGER NPI , NPM , NPMP J , NP 1 

REAL Cl C2 ,C3 , DELI , DELXM, DELXMM, DELXl , DELX2 ,DELYN,DBLYNM, DELYl 
Seal Sely2;dbli,del2,diagi,diagim,diag1,diag2,one,sdiagi,sdiag2 
REAL T,ZER0,ZXYMNS,ZXY1NS 

CONSTRUCT CONSTANTS 

data zero, one/o. OEO, 1 .oeo/ 

DEFINE SOME VARIABLES 


MMl - 
MPl = 
NM1 * 
NPl * 
NPM * 
IERR 


M - 1 
M + 1 
N - 1 
N + 1 
N + M 
= 2 


C 

c 


OBTAIN Y-PARTIAL DERIVATIVES ALONG Y = Y(l) 

DELYl = DELY „ ,, 

IF (IOPTY .EQ. 0) DELYl = Y(2) - Y(l) 

IF (DELYl .LE. ZERO) GO TO 9000 
IF ( SIGMAY*DELY 1 . LE . 650.0) GO TO 5 
IERR » 3 
GO TO 9000 

5 IF (ISLPSW(3) .EQ. 2) GO TO 100 
IF ( ISLPSW( 3 ) .EQ. 1) GO TO 20 

THIS CODE IS FOR A USER-DEFINED EDGE CONDITION 



c 

c 

c 


c 

c 

c 


DO 10 1 « lfM 

zf(l,l,l) - zYi(t) 

10 CONTINUE . 4 

IF ( ISLPSW( 4 ) . NE . 1) GO TO 100 
IF (IOPTY .EQ. 0) GO TO 140 

THIS CODE IS FOR A COMPUTER-DEFINED EDGE CONDITION 

20 « .GT. 2, DELV2 - Y(3. - »d. 

IF (DELY2 .LE. DELYl) GO TO 9000 

CALL QXZ060 ( DELYl , DELY2 , SIGH AY ,Cl f C2,C3 r N) 

BRANCH SINCB ONLY THE Cl’S ABE NEEDED 

IF (ISLPSWp) .EQ. 0) GO TO 130 

DO 30 I - l r M 

ZP(I,1,1) » Cl*Z(I,l) + C2*Z(I,2) 

30 CONTINUE 

IF (N .EQ. 2) GO TO 100 
DO 40 I = 1,M 


ZP(I, 1, 1) 
40 CONTINUE 


ZP(I,1,1) + C3*Z ( 1 , 3 ) 


C 

C 

C 


OBTAIN Y-PARTIAL DERIVATIVES ALONG Y = Y(N) 

100 IF ( ISLESW( 4 ) .EQ. 2) GO TO 200 
IF ( I SLPSW ( 4 ) .EQ. 1) GO TO 120 


C THIS CODE IS FOR A USER-DEFINED EDGE CONDITION 
DO 110 I - 1,M 
TEMP(N+I) - ZYN(I) 

110 CONTINUE 
GO TO 200 

C THIS CODE IS FOR A COMPUTER-DEFINED EDGE CONDITION 
120 IF (IOPTY .EQ. 0) GO TO 140 

C THIS CODE IS FOR THE Y-EQUALLY-SPACED CASE 
C 

130 Cl - -cl 
C2 - -C2 

IF (N .GT. 2 ) C3 = -C3 
GO TO 150 


C 

C 

C 


C 

C 


THIS CODE IS FOR THE NON-Y-EQUALLY-SPACED CASE 

140 DELYN - Y ( N ) - Y(NMl) 

IF (DELYN .LE. ZERO) GO TO 9000 
IF ( SIGMAY*DELYN .LE. 650.0) GO TO 145 
IERR = 3 
GO TO 9000 

145 DELYNM - DELYN + DELYN 

IF (N .GT. 2) DELYNM = Y(N) - Y(N-2) 

IF (DELYNM .LB. DELYN) GO TO 9000 

CALL QXZ060 ( -DELYN , -DELYNM , SIGMAY ,Cl ,C2 ,C3 ,N) 

150 TEMP(N+I) - Cl*Z( I ,N) + C 2 * Z ( I , NM 1 ) 

160 CONTINUE 

IF (N -EQ. 2) GO TO 200 
DO 170 I - l f M 
NPI * N + I 

TEMP(NPI) - TEMP(NPI) + C3*Z(lrN-2) 

170 CONTINUE 

OBTAIN X-PARTIAL DERIVATIVES ALONG X = X(l) 



c 


200 bELXl - DfeLX 

ir (loETX .EQ. 0) bELXl = X{2) - X ( 1 ) 

IF (DELXl .LB. ZERO) GO TO 9000 
IF ( S IGM AX * DE LX 1 .LE. 650.0) GO TO 205 
IERR « 3 
GO TO 9000 

205 IF (ISLPSW(l) .EQ. 2) GO TO 400 

IF (ISLPSW(l) .EQ. 1) GO TO 220 

C 

C THIS CODE IS FOR A USER-DEFINED EDGE CONDITION 
C 

DO 210 J - 1 , N 
ZP(1,J,2) = ZXl(J) 

2 10 CONTINUE 

IF 1 ISLPSW( 5 ) .EQ. 1) GO TO 220 

IF ( ISLPSW( 7 ) .EQ. 1) GO TO 220 

IF (IOPTX .EQ. 0) GO TO 310 
IF (ISLPSW(2) .EQ. 1) GO TO 220 

IF (ISLPSW(6) .EQ. 1) GO TO 220 

IF ( I SLPSW ( 8 ) .EQ. 0) GO TO 310 

C 

C THIS CODE IS FOR A COMPUTER-DEFINED EDGE CONDITION 
C 

220 DELX2 = DELXl + DELXl 

IF ( IOPTX .EQ. 0 .AND. M .GT. 2) DELX2 => X(3) - X(l) 

IF (DELX2 .LE. DELXl) GO TO 9000 

CALL QXZ060( DELXl , DELX2 , SIGMAX , Cl , C2 ,C3 ,M) 

C 

C BRANCH SINCE ONLY THE Cl'S ARE NEEDED 

C 

IF (ISLPSW(l) .EQ. 0) GO TO 300 
DO 230 J - 1,N 

ZP( 1 , J, 2 ) - Cl*Z(l,J) + C2*Z( 2, J) 

230 CONTINUE 

IF (M .EQ. 2) GO TO 300 
DO 240 J = 1,N 

ZP( 1, J r 2 ) - ZP( 1 , J, 2 ) + C3*Z { 3 , J ) 

240 CONTINUE 
C 

C OBTAIN X-Y-PARTIAL DERIVATIVE AT ( X ( 1 ) , Y ( 1 ) ) 

C 

300 IF ( ISLPSW( 5 ) .EQ. 1) GO TO 320 
C 

C THIS CODE IS FOR A USER-DEFINED CORNER CONDITION 
C 

310 ZP( 1, 1,3) - ZXYll 
GO TO 330 


C 

C 

C 


C 

C 

C 

C 

c 


THIS CODE IS FOR A COMPUTER-DEFINED CORNER CONDITION 


320 


ZP( 1 , 1 , 3 ) 
IF (M .GT. 


* Cl*ZP( 1,1,1) 
2) ZP( 1,1,3) 


+ C2*ZP( 2, 1 , 1 ) 

- ZP( 1 , 1 , 3 ) + C3*ZP( 3,1,1) 


OBTAIN X-Y-PARTIAL DERIVATIVE AT (X(1),Y(N)) 


330 IF ( ISLPSW( 7 ) .EQ. 1) GO TO 340 

THIS CODE IS FOR A USER-DEFINED CORNER CONDITION 
ZXY1NS = ZXY1N 
GO TO 500 


C 

C THIS CODE IS FOR A COMPUTER-DEFINED CORNER CONDITION 

C 

340 ZXYlNS - Cl*TEMP(N+l ) + C2*TEMP(N+2) 

IF (M .GT. 2) ZXYlNS = ZXYlNS + C3*TEMP(N+3) 



GO TO 500 

C THIS CObE IS PoR A NATURAL EDGE ALONG X - X(l) 
C 

400 DO 410 J=1,N 

ZP(1,J,2) - ZERO 
410 CONTINUE 

ZP(1,1,3) - ZERO 
ZXYlNS - ZERO 


C 

C 

C 


C 

C 

C 


C 

C 

C 

C 

C 

C 


OBTAIN X-PARTIAL DERIVATIVE ALONG X - X(M) 

500 IE ( ISLPSW( 2 ) .EQ. 2) GO TO 700 

IP ( ISLPSW( 2 ) .EQ. 1) GO TO 520 

THIS CODE IS POR A USER-DEPlNED EDGE CONDITION 

DO 5 10 J = 1,N 
TEMP(NPM+J) - ZXM ( J ) 

510 CONTINUE ^ 

IF (ISLPSW(6) .EQ. 1) GO TO 520 

IF (ISLPSW(8) .EQ. 0) GO TO 610 

THIS CODE IS FOR A COMPUTER-DEFINED EDGE CONDITION 

520 IF (loPTX .EQ. 0) GO TO 530 

THIS CODE IS POR THE X-EQUALLY-SPACED CASE 


Cl - -cl 
C2 - -C2 

IF (M .GT. 2) C3 
GO TO 540 


-C3 


C 

C 

C 


c 

c 

c 


THIS CODE IS POR THE NON-X-EQUALLY-SPACED CASE 

530 DELXrt - X(M) - X(MM1) 

IF (DELXM .LE. ZERO) GO TO 9000 
IF (SIGMAX*DELXM .LE. 650.0) GO TO 535 
1ERR = 3 
GO TO 9000 

535 DELXMM - DELXM + DELXM u 

IF (M .GT. 2) DELXMM - X(M) - X(M-2) 

IF (DELXMM .LE. DELXM) GO TO 9000 

CALL QXZ060 ( -DELXM, -DELXMM, SIGMAX ,Cl ,C2 ,C3 ,M) 

BRANCH SINCE ONLY THE Cl'S ARE NEEDED 

540 IF (ISLPSW(2) .EQ. 0) GO TO 600 

DO 550 J « 1,N 

TEMP(NPM+J ) - Cl*Z(M,J) + C2*Z(MM1,J) 

550 CONTINUE 

IF (M .EQ. 2) GO TO 600 
DO 560 J * 1 , N 
NPMPJ ” NPM + U 

TEMP(NPMPJ) - TEMP ( NPMPJ ) + C3*Z(M-2,J) 

560 CONTINUE 


C OBTAIN X-Y-PARTIAL DERIVATIVE AT (X(M),Y(1)) 

C 600 IF ( ISLPSW( 6 ) .EQ. 1) GO TO 620 
C THIS CODE IS FOR A USER-DEFINED CORNER CONDITION 

C 610 ZP(M,1,3) - ZXYM1 

GO TO 630 



non non non nnn 


tills CODE is toll A COMPUTER-DEFINED CORNER CONDITION 

cm 7 t>(H | it b cl*ZP(M, 1,1) + C2*ZP(MMi , 1,1) 

620 ,, JP ,i,S,3) . SP(M,1,3| * C3.SP(M-2,1,1) 

OBTAIN X-Y-PARTlAL DERIVATIVE At (X(M),Y(N)) 

630 It (1SLPSW(8) .EQ. 1) GO TO 640 

THIS CODE IS FOR A USER-DEFINED CORNER CONDITION 

ZXYMNS = ZXYMN 
GO TO 800 

THIS CODE IS FOR A COMPUTER-DEFINED CORNER CONDITION 

c4n 7XYMNS » Cl*TEMP(NPM) + C2*TEMP(NPM-1 ) 

640 IF (M.GT. 2) ZXYMNS ■ ZXYMNS + C3*TEMP(NPM-2 ) 

GO TO 800 

; this CODE IS FOR A NATURAL EDGE ALONG X ** X(M) 

"I 

700 DO 710 >1=1 , N 

TEMP(NPM+J) = ZERO 
710 CONTINUE 

ZP(M, 1,3) = ZERO 
ZXYMNS = ZERO 

r CPT UP RIGHT HAND SIDES AND TRIDIAGONAL SYSTEM FOR Y-GRID 
c PERFORM FORWARD ELIMINATION ON THE FIRST COLUMN 

c 

000 DELl = DELYl 

DELI = ONE/DELl 

DO 810 I = 1 ,M - f »ni 

ZP( 1 , 2 , 1 ) = DELIMZ(lr2) - Z(I,1)) 

7P( 1 2.3) = DELI * ( ZP( 1 , 2 , 2 ) - ZP(1,1,2)) 

ZP( m! 2 , 3 ) = DELI * ( TEMP ( NPM+2 ) - TEMPCNPMPl ) ) 

CALL QXZ061 ( D1AG 1 , SDIAGl , S1GMAY ,DEL1 ) 

DIAGI = ONE/DIAG1 
IF (ISLPSW(3) .EQ. 2) GO TO 830 

C THIS CODE IS FOR A NON-NATURAL EDGE ALONG Y = Y(l) 

ZP< 1 , 1 , 3 ) = DIAGI*(ZP( 1,2,3) - JojA'l',!! 

ZP(M, 1,3) = DIAGI*(ZP(M,2,3) - ZP(M,1,3)) 

TEMP ( 1 ) = DIAGI *SDIAGl 

DO 820 I = 1 , M , ... 

ZP(I,1,1) = DIAGI * ( ZP ( 1 , 2 , 1 ) - ZPd,l,D) 

820 CONTINUE 
GO TO 850 

C THIS CODE IS FOR A NATURAL EDGE ALONG Y = Y(l) 

C 

830 DO 840 1=1, M 

ZP(I,1,1) = ZERO 
840 CONTINUE 

ZP ( 1 , 1 , 3 ) = ZERO 
ZP(M, 1,3) = ZERO 
TEMP(l) = ZERO 
850 IF (N .EQ. 2) GO TO 1000 

C THESE TERMS ARE DEFINED IN THE Y- EQUALLY -SPACED CASE 


DEL2 = DELl 



non nrao nnoo 


DIAG2 - blAGi 
SblAG2 * SblAGl 

THIS IS tMfi FoRWARb ELIMINATION Loot 
LET THE WILD RUMPUS BEGIN 

DO 930 J = 2 ,NMl 
JM1 = J - 1 
JP1 *» J 4- 1 
NPMPJ « NPM + J 
IF ( IOPTY .NE. 0) GO TO 900 

THIS CODE IS FOR THE NON-Y- EQUALLY- SPACED CASE 

DEL2 - Y(JP1) - Y(J) 

IF (DEL2 .LE. ZERO) GO TO 9000 
IF (SIGMAY*DEL2 .LE. 650.0) GO TO 890 
IERR = 3 
GO TO 9000 
890 DELI = ONE/DEL2 

CALL QXZ06 1 ( DIAG2 , SDI AG2 , SIGMAy , DEL2 ) 

CONTINUE THE FORWARD ELIMINATION LOOP 

900 DO 910 I = 1 ,M 

ZP(I,JPl,l) - DELI * ( Z ( I , JPl ) - Z(I,J)) 

910 CONTINUE 

ZP( 1 , JPl , 3 ) - bELIMZP( - ZP( 1 , J, 2 ) ) 

ZPlM, JPl , 3 j - DELI* (TEMP(NPMPJ+1 ) - TEMP(NPMPJ)) 
DIAGIN - ONE/ ( DI AG 1 + DIAG2 - SDIAGl*TEMP( JM1 ) ) 

DO 920 I - 1 f M 

ZP(I,J,1) - DIAGIN* (ZP(I»JPlfl) - ZP(I,J,1) - 

* SDIAGl*ZP( I, JM1 , 1 ) ) 

920 CONTINUE 

ZP(l,J f 3) - DIAGIN*(ZP( l,JPl,3) - ZP(1,J,3) - 

* SDIAG1*ZP( 1, JM1, 3) ) 

ZP(M, J,3) - DIAGIN* (ZP(M,JPl f 3) - ZP(M,J,3) - 

* SDIAGl*ZP(M, JMl r 3) ) 

TEMP(J) = DIAGIN* SDI AG2 
DIAG1 = DIAG2 
SDI AG 1 = SDIAG2 

930 CONTINUE 
C 

C END OF THE FORWARD ELIMINATION LOOP 
C NOW PERFORM FORWARD ELIMINATION ON THE LAST COLUMN 

C 

1000 IF ( ISLPSW( 4 ) .EQ. 2) GO TO 1020 

C 

C THIS CODE IS FOR A NON-NATURAL EDGE ALONG Y - Y(N) 

c 

DIAGIN = ONE/(DlAGl - SDIAGl *TEMP( NM1 ) ) 

ZP( 1 ,N, 3 ) = DIAGIN* ( ZXY INS - ZP(1,N,3) - 

* SDIAGl *ZP ( 1 , NM1 , 3 ) ) 

TEMP(N) - DIAGIN* (ZXYMNS - ZP(M,N,3) - 

* SDIAGl *ZP(M f NMl , 3 ) ) 

DO 1010 I = 1,M 

ZP{ I,N, 1 ) - DIAGIN* (TEMP(N+I ) - ZP(I,N,1) - 

* SDIAGl *ZP ( I , NM1 , 1 ) ) 

1010 CONTINUE 

GO TO 1100 

C 

C THIS CODE IS FOR A NATURAL EDGE ALONG Y = Y(N) 

C 

1020 DO 1030 1=1, M 

ZP( I ,N, 1 ) = ZERO 
1030 CONTINUE 



n n n o n 


ZP(l,N r 3) - ZERO 
TEMP(N) « ZERO 

PERFORM BACK SUBSTITUTION 

1100 DO 1120 J • 2 , N 
JBAK - NPl - J 
JBAKPl = JBAK + 1 
T - TEMP (JBAK) 

ZP( I ) JBAK , 1 ) - ZP ( I , JBAK , 1 ) - T*ZP( I , JBAKPl, 1 ) 

mo continue _ ZP(1/JBAKf3 , - T*zf(i, jbakpi,3) 

TEMP(JBAK) = ZP ( M, JBAK , 3 ) - T*TEMP( JBAKPl ) 

1120 CONTINUE 

SFT UP RIGHT HAND SIDES AND TRIDIAGONAL SYSTEM FOR X-GRID 

S»lSS a NATURAL EDGE 


DELI - DELXl 
DELI = ONE/DELl 
DO 1200 J ■ 1»N 
ZP(2,J,2) ** BELI*(Z(2,J)^ 


Z(1,J) 


ZP( 2 , J , 3 j - DELI * ( ZP ( 2 , J , 1 ) - ZP 
1200 CONTINUE 


!U 


m 


CALL QXZ061 (DlAGl r SDIAGl , SIGMAX ,DELl ) 

DIAGI - ONE/DIAGl 
IF (ISLPSW(l) . EQ. 2) GO TO 1220 

ZP( 1^ J^2 > ) ** DIAGI * ( ZP ( 2 , J , 2 ) - ZP(1,J,2)| 

ZP(l|j,3) « DIAGI * ( ZP( 2 , J , 3 ) - ZP(1,J.3)) 

1210 CONTINUE 

TFMPf NPl 1 = dxAgx*sdiagi 
1220 i?(iSSi*(l) .EQ. 2j TEMP(NPl) = ZERO 
IF (M .EQ. 2) GO TO 1400 

THESE TERMS ARE DEFINED IN THE X- EQUALLY -SPACED CASE 

DEL2 = DELl 
DIAG2 = DIAGI 
SDIAG2 = SDIAGl 


C 

C 

C 


C 

C 

C 


C 

C 

c 


c 

c 

c 


THIS IS THE FORWARD ELIMINATION LOOP 

DO 1330 I = 2 ,MMl 
IM1 - I - 1 
IPl =1+1 
HPI “ N + X 

IF ( lOPTX .NE. 0) GO TO 1300 
THIS CODE IS EOR THE NON-X-EQUALLY-SPACED CASE 

DEL2 - X(IPl) - X(I) 

IF (DEL2 .LE. ZERO) GO TO 9000 
IF (SIGMAX*DEL2 .LE. 650.0) GO TO 1290 

I ERR « 3 
GO TO 9000 

12,0 CALL QXZ06 l^DIAG2 , SDIAG2 , SIGMAX,DEL2 ) 
CONTINUE THE FORWARD ELIMINATION LOOP 

1300 zp(ipi?j?2) *='DELI*(Z(IP1,J) - ZU,JJ) 

ZP(lPl»J»3) « DELI * ( ZP ( IP 1 , J t 1 ) “ ZP(I,J,1)) 



non n n o a 


1310 CONTINUE 
DlAGIN - 
DO 1320 3 
ZP(I,J,2) 

* 

ZP( 1 , J r 3 ) 


ONE/(DlAGl + blAG2 - 

- 1,N 

- DIAGlN*(ZP(IPl,«J, 2) 
SDlAGl*ZP(lMl f J f 2)) 

* D I AG IN* (ZP(IPlfJf3) 
SDlAGl*ZP( IM1 , J, 3) ) 


1320 - DIRG1M*SDIAG2 

DIAGl = DIAG2 
SDIAGl = SDIAG2 
1330 CONTINUE 


SDIAG1*TEMP(NPI - 1)) 

- ZP(I,J,2) - 

- ZP(I,J,3) - 


END OF THE FORWARD 
NOW PERFORM FORWARD 


ELIMINATION LOOP 

elimination on the 


LAST COLUMN 


1400 IF ( 1SLPSW( 2 ) .EQ. 2) GO TO 1420 

THIS CODE IS FOR A NON-NATURAL EDGE ALONG X = X(M) 


DlAGIN = ONE/ (DIAGl - SDIAGl*TEMP(NPM - D) 


DO 1410 J 

ZP(M, J, 2 ) 

ZP(M, J,3) 


ZP(M, J, 2 ) 


1,N 

DlAGIN* (TEMP(NPM+J) 
SDIAGl*ZP(MMl, J,2) ) 

DlAGIN* (TEMF(J) - ZP(M, J , 3 ) - 
SDIAGl *ZP( MM 1 , J, 3) ) 


1410 CONTINUE 

GO TO 1500 


C THIS CODE IS FOR A NATURAL EDGE ALONG X 
C 

1420 DO 1430 J=1,N 

ZP(M,J,2) - ZERO 
ZP(M,J,3) - ZERO 
1430 CONTINUE 


C 

C 

C 


PERFORM BACK SUBSTITUTION 


1500 DO 1520 I = 2,M 
IBAK = MPl - I 
IBAKPl ” IBAK + 1 
T = TEMP (N+ IBAK) 

DO 1510 3 = IfN 
ZPC IBAK, J t 2 ) - ZP ( IBAK , 
ZP( IBAK f J , 3 ) = ZP ( IBAK , 
1510 CONTINUE 
1520 CONTINUE 
IERR = 0 
9000 RETURN 
END 

GO TO 1500 


3 , 

3 , 


2 ) 

3) 


T*ZP(IBAKPl, J,2) 
T*ZP( IBAKPl » Jr 3) 


C THIS CODE IS FOR A NATURAL EDGE ALONG X = X(M) 

c 

1420 DO 1430 J=1,N 

ZP(M, J, 2 ) = ZERO 
ZP(M,J,3) = ZERO 
1430 CONTINUE 

C PERFORM BACK SUBSTITUTION 

c „ 

1500 DO 1520 I = 2 , M 

IBAK » MPl - I 
IBAKPl = IBAK + 1 
T = TEMP(N+IBAK) 

DO 1510 J - 1 • N 



